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ABSTRACT 

Simple assumptions made regarding electron thermodynamics often limit the extent to which 
general relativistic magnetohydrodynamic (GRMHD) simulations can be applied to obser¬ 
vations of low-luminosity accreting black holes. We present, implement, and test a model 
that self-consistently evolves an entropy equation for the electrons and takes into account the 
effects of spatially varying electron heating and relativistic anisotropic thermal conduction 
along magnetic field lines. We neglect the back-reaction of electron pressure on the dynamics 
of the accretion flow. Our model is appropriate for systems accreting at <sc 10“^ of the Edding¬ 
ton accretion rate, so radiative cooling by electrons can be neglected. It can be extended to 
higher accretion rates in the future by including electron cooling and proton-electron Coulomb 
collisions. We present a suite of tests showing that our method recovers the correct solution 
for electron heating under a range of circumstances, including strong shocks and driven tur¬ 
bulence. Our initial applications to axisymmetric simulations of accreting black holes show 
that (1) physically-motivated electron heating rates that depend on the local magnetic field 
strength yield electron temperature distributions significantly different from the constant elec¬ 
tron to proton temperature ratios assumed in previous work, with higher electron temperatures 
concentrated in the coronal region between the disc and the jet; (2) electron thermal conduc¬ 
tion significantly modifies the electron temperature in the inner regions of black hole accretion 
flows if the effective electron mean free path is larger than the local scale-height of the disc 
(at least for the initial conditions and magnetic field configurations we study). The methods 
developed in this work are important for producing more realistic predictions for the emis¬ 
sion from accreting black holes such as Sagittarius A* and M87; these applications will be 
explored in future work. 
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1 INTRODUCTION 

A wide variety of low luminosity accreting black holes are cur¬ 
rently interpreted in the context of a Radiatively Inefficient Ac¬ 
cretion Flow (RIAF) model that describes a geometrically thick, 
optically thin disc with a low accretion rate and luminosity. In par¬ 
ticular, this is true of the black hole at the center of our galaxy, 
Sagittarius A* (Narayan et al. 1998), the black hole at the center 
of Messier 87 (Reynolds et al. 1996), and other low luminosity Ac¬ 
tive Galactic Nuclei (AGN), as well as a number of X-ray binary 
systems (see Remillard & McClintock 2006 for a review). The gas 
densities in these systems are low enough that the time scale for 
electron-ion collisions is much longer than the time scale for accre¬ 
tion to occur, so a one-temperature model of the gas is no longer 
valid (as originally recognised by Shapiro, Lightman & Eardley 
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1976, Ichimaru 1977, and Rees, Begelman, Blandford, & Phinney 
1982). Instead, a better approximation is to treat the electrons and 
ions as two different fluids, each with its own temperature. 

Calculating the emission from accreting plasma requires pre¬ 
dicting the electron distribution function close to the black hole. To 
date, time dependent numerical models of RIAFs that attempt to 
directly connect to observations often assume a Maxwellian distri¬ 
bution with a constant electron to proton temperature ratio, Tg/Tp, 
and take the results of GRMHD simulations as the solution for the 
total gas temperature, Tg = Tp + Tg (Dibi et al. 2012; Drappeau et al. 
2013; Moscibrodzka et al. 2009). This neglects, however, several 
physical processes that have different effects on the electron and 
proton thermodynamics and that are currently only included in one¬ 
dimensional semi-analytic models. Such effects include electron 
thermal conduction (e.g., Johnson & Quataert 2007), electron cool¬ 
ing (e.g., Narayan & Yi 1995), and non-thermal particle accelera¬ 
tion and emission (e.g.. Yuan, Quataert & Narayan 2003). To date. 
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extensions of the simple TpjT,, = const, prescription have been lim¬ 
ited to post-processing models that do not self-consistently evolve 
the electron thermodynamics over time. Examples include the pre¬ 
scription of Moscibrodzka et al. (2014) which takes TpIT^ = const, 
in the disc proper but sets Tg = const, in the jet outflow region, 
as well as the model of Shcherbakov, Penna & McKinney (2012), 
who solve a 1-D radial equation for Tp - at a single time-slice in 
the midplane to obtain a functional relationship TpITg = f{Tg) that 
is then applied to the rest of the simulation. To enable a more ro¬ 
bust connection between observations of accreting black holes and 
numerical models of black hole accretion, it is critical to extend the 
detailed thermodynamic treatment of electrons used in ID calcula¬ 
tions to multi-dimensional models. This is the goal of the current 
paper. In particular, we describe numerical methods for separately 
evolving an electron energy equation in GRMHD simulations. We 
focus on including heating and anisotropic thermal conduction in 
these models. Future work will include electron radiative cooling 
and Coulomb collisions between electrons and protons. 

In a turbulent, magnetised plasma, electrons and ions are 
heated at different rates depending on the local plasma conditions 
(e.g., Quataert & Gruzinov 1999; Cranmer et al. 2009; Howes 
2010; Sironi 2015). Furthermore, since the electron-to-proton mass 
ratio is small, electrons will both conduct and radiate their heat 
much more efficiently than the ions. The combination of these ef¬ 
fects leads to the expectation that, in general, Tg < Tp. In the 
present paper, we thus neglect the effect of the electron thermody¬ 
namics on the overall dynamics of the accretion flow. This allows 
us to treat the simulation results as a fixed background solution on 
top of which we independently evolve the electrons. Even if we find 
that Tg ~ Tp in some regions of the disc, this treatment may still 
be a reasonable first approximation given the uncertainties in the 
electron physics. 

The neglect of electron cooling in the present paper is reason¬ 
able for systems accreting at < 10“^ of the Eddington rate, M^u, so 
that the synchrotron cooling time is much longer than the accretion 
time (Mahadevan & Quataert 1997). In particular, this likely in¬ 
cludes Sagittarius A* in the galactic center. The application of our 
methodology to Sgr A* is particularly important given the wealth 
of multi-wavelength data (e.g., Serabyn et al. 1997, Zhao et al. 
2003, Genzel et al. 2003, Baganoff et al. 2003, Barriere et al. 2014) 
and current and forthcoming spatially resolved observations with 
the Event Horizon Telescope (Doeleman et al. 2008) and Gravity 
(Gillessen et al. 2010). 

The goal of this paper is to present our formalism and method¬ 
ology for evolving the electron thermodynamics and to apply the 
results to 2D (axisymmetric) GRMHD simulations of an accreting 
black hole. We show the range of possible electron temperature dis¬ 
tributions in the inner region of the disc, which directly impacts the 
predicted emission. Future work will explore the impact that these 
results have on the emission, spectra, and images of Sagittarius A*. 

The remainder of this paper is organised as follows. §2 de¬ 
scribes our theoretical model of electron heating and anisotropic 
electron conduction while §3 describes the numerical implementa¬ 
tion of this model. §4 contains tests of the numerical implementa¬ 
tion, §5 applies the model to a 2D simulation of an accretion disc 
around a rotating black hole, and §6 discusses the implications of 
this application and concludes. Boltzmann’s constant, fc*, and pro¬ 
ton mass, nip, are taken to be 1 throughout. We use cgs units, with 
Lorentz-Heaviside units for the magnetic field (e.g., magnetic pres¬ 
sure is E12), and a metric signature of (—l-l-l-). We also assume 
that the gas is mostly hydrogen and ideal. Since we also assume 


Hg X Up = n, then p = rUgUg + nipHp x nipn = n (setting nip = 1), so 
we use p and n interchangeably. 


2 ELECTRON THERMODYNAMICS 

The accreting plasmas of interest are sufficiently low density that 
the electron-proton Coulomb collision time is much longer than the 
dynamical time and so a two-temperature structure can develop, 
with the protons and electrons having different temperatures (and, 
indeed, different distribution functions). Moreover, at the low ac¬ 
cretion rates where radiative cooling can be neglected, the electron- 
electron and proton-proton Coulomb collision times are also much 
longer than the dynamical time (Mahadevan & Quataert 1997). 
However, the plasma densities are high enough that the plasma 
is nearly charge-neutral and so we assume that Ug x np. We fur¬ 
ther assume that the electron flow velocity is the same as that of 
the protons.* This need not strictly be true (e.g., in the solar wind 
the relative velocities of particle species can be of order the Alfven 
speed; e.g. Bourouaine et al. 2013 ), but is a reasonable first approx¬ 
imation. A similar approach is often used in modeling the global 
dynamics of the low-collisionality solar wind (e.g., Chandran et al. 
2011 ). 

Under these assumptions, the key difference in the electron 
and proton physics lies in their different thermodynamics: the pro¬ 
tons and electrons have very different heating and cooling pro¬ 
cesses that need to be separately accounted for. Formally, because 
of the low collisionality conditions we should separately solve the 
electron and proton Vlasov equations. This is computationally ex¬ 
tremely challenging, however, particularly in the global geometry 
required to predict the emission from accreting plasmas (even lo¬ 
cal shearing box calculations using the particle-in-cell technique 
to solve the Vlasov equation require an unphysical electron-proton 
mass ratio, thus making it difficult to reliably model the electron 
thermodynamics; e.g., Riquelme et al. 2012). As a result, we as¬ 
sume a fluid model in this paper. Our fluid approximation corre¬ 
sponds to taking moments of the Vlasov equation and applying 
closures on higher moments of the distribution function. As we 
shall describe, our closure corresponds to specific models for the 
conductive heat flux, the viscous momentum flux, and the turbulent 
heating rate of each particle species. 

Our basic model is thus to take a single-fluid GRMHD so¬ 
lution (e.g., Komissarov 1999; Gammie, McKinney & Toth 2003; 
De Villiers & Hawley 2003) as an accurate description of the total 
fluid (composed of both the electron and proton gas) dynamics and 
thus the accretion flow density, magnetic field strength, and velocity 
field. We evolve the electrons as a second fluid on top of this back¬ 
ground solution. The GRMHD solution may itself include viscosity 
and conduction as in Chandra et al. (2015). Our assumption that the 
electrons do not back react on the flow dynamics is formally valid 
in the limit that Tg « Tp, but should be a reasonable approximation 
so long as Tg < Tp in regions of large plasma j8 > 1, i.e., where gas 
pressure forces are dynamically important. One advantage of not 
coupling the electron pressure to the GRMHD solution is that we 
can run multiple electron models in one simulation, allowing us to 
explore systematic uncertainties with a minimum of computational 
time. 

In this paper, we focus on implementing electron heating and 

* More precisely, as in standard MHD, the relative velocity between elec¬ 
trons and protons required to produce currents that can maintain magnetic 
fields near /? ~ 1 is orders of magnitude less than the mean sound speed. 
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anisotropic conduction. Coulomb collisions are straightforward to 
include but are negligible for the low accretion rates at which elec¬ 
tron cooling can be neglected. In future work, electron cooling will 
be self-consistently incorporated building on the BHlight code de¬ 
veloped by Ryan, Dolence & Gammie (2015). 


2.1 Basic Model 

The stress-energy tensors for the electron and proton fluids in our 
model take the form: 

rf ^{p, + u, + P ,)« -H + + + 

rr = (pp + u, + P,) + T^;, 

wherep*,, uj,, and P*. are the fluid frame density, internal energy, and 
pressure, respectively, is the fluid four-velocity in the coordinate 
frame, is a general stress tensor that accounts for viscous effects, 
and is the heat flux carried by the electrons. The subscript k 
denotes p or e (and will also denote the total gas quantities labeled 
by g below). We leave as a general tensor that will be model- 
specific. For each species, ignoring electron-electron, electron-ion, 
and ion-ion collisions, one can take the zeroth and first moment of 
the Vlasov equation to show that 

Vp (p*<) = 0 (2) 

and 

(3) 

Vpr;’' = 

where is the electromagnetic field tensor. In ideal, single-fluid 
GRMHD in the absence of shocks, the conservation of entropy 
equation, pTgU>^dijSg = 0, where Sg is the entropy per particle, 
follows directly from the conservation of particle number and the 
stress-energy (see page 563 in Misner, Thome & Wheeler 1973). 
To derive entropy equations for the electron and proton fluid used 
in our model, we perform the same series of manipulations; namely, 
contracting both equations (3) with u'' (the total fluid velocity, 
which we take to be =» » u^) and invoking equation (2), which 

give us: 

pTXdi^Sp = Qp, (4) 

and 

pTXdpS, = Q,-Vpq>l-apif^, (5) 

where we have defined the heating rate per unit volume for each 
species as a sum of viscous and Ohmic resistance terms, Qg = 
UvVpT^'' -I- enUl,u''Fp and Qp = UyVpi^'' - eni/pU''Fp, and where 
fl'' = is the four-acceleration, which accounts for gravita¬ 

tional redshifting of the temperature by the metric. We can write the 
heating rates in terms of the electric field four-vector, = UyF*''', 


and the four-currents, 7^ = —nedl, Jp = neiTp as^: 

Qy = UyVpF: + J^ep 

Qp = UyVpF^p'' + FpBp. 

The intuitive understanding of the Ohmic heating terms on the 
right-hand side of equation (6) is that they are ~ Jh ■ E evaluated 
in the rest frame of the total fluid. To derive the entropy equation 
for the total fluid, we first define several total fluid variables as a 
sum of electron and proton terms: p = Pp + pc x n,Ug = Up + u^, 
Pg = Pp + Pc, Tg = Tp + Tc, Ji‘ = Tp + Jl = en{i/p - lil) and = 
FT -H FT, denoting total gas mass density, internal energy, pressure, 
temperature, current, and viscous stress. Then, using the thermody¬ 
namic identity, pT^u^dpSk = iFdpUjc - (uk + PQ u^dp log(p), we find 
that the entropy per particle of the total gas, Sg, satisfies the relation 
pTguTdpSg = pTpuTdpSp + pTcuTdpSc, resulting in 

pTgiFdpSg = Q- Slpcf^ - ap(f^, (7) 

with the total heating rate per unit volume: 

Q=Q„ + Qe = UyVpFT + Eep. ( 8 ) 

In practice, we use the electron entropy equation (5) to evolve 
the electron thermodynamics. To determine the overall dynamics of 
the electron -r proton gas, we use Maxwell’s equations in addition 
to a standard, single-fluid GRMHD evolution representing the total 
gas. The equations for the latter are obtained by separately sum¬ 
ming the electron and proton parts of equation (2) and equation (3), 
resulting in a mass conservation equation. 


Vp (piF) = 0, 

(9) 

and an energy-momentum equation. 



(10) 

with the total gas stress-energy tensor. 


TT^{p + Ug + Pg)iFu'’ + Pgf\ 

(11) 


and the electromagnetic stress energy tensor^, = F^°‘F^ — 
g^''FafiF‘‘^ ! A. We have used the identity '^TF'^Tm ~ equa¬ 

tion (10), the assumption that i4 ~ ufl, ~ iF in equation (11), and 
the charge-neutrality assumption iic = rip = n throughout. Further¬ 
more, we have dropped the electron thermal conduction terms in 
the evolution of the total gas properties, though we keep them in the 
evolution of the electron entropy (equation 5). This is consistent if 
Tc < Tp since electron conduction will then affect the electron ther¬ 
modynamics but not the overall stress-energy of the fluid. Finally, 
we take T^^ to be given by the ideal MHD limit (i.e., F —> 0): 

= + ( 12 ) 

where is the magnetic field four-vector defined 

^ We have kept the subscripts e and p for the four-currents (and thus 
four-velocities) in equation (6) because the details of the Ohmic heating 
depend on the small but non-zero velocity difference between the proton 
and electron fluids (or equivalently the velocity difference between the elec¬ 
tron/proton fluid and the total fluid). An explicit expression for these terms 
would require a detailed kinetic theory calculation beyond the scope of the 
present work (i.e., some form of “generalized Ohm’s Law,” as in, e.g., Koide 
2010). In our model, as described in the text, numerical resistivity provides 
the Ohmic heating that is then distributed to electrons and protons according 
to a closure model obtained from previous work in kinetic theory. 

^ Here we have chosen to absorb a factor of (4;r)“^^^ into the definition of 
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in terms of the Levi-Civita tensor, and E = is twice the 
magnetic pressure. With these assumptions, equations (9), (10), and 
Maxwell’s equations are simply the standard single-fluid equations 
of ideal GRMHD except with an explicit viscosity tensor. In stan¬ 
dard conservative GRMHD codes (including the one used in this 
work), this viscosity tensor is not included explicitly hut implic¬ 
itly generated numerically hy the Riemann solver. Furthermore, the 
Riemann solver also introduces a finite numerical resistivity into 
Maxwell’s equations, allowing for a nonzero (and thus nonzero 
Ohmic heating). For further discussion of these points, see § 3.1. 

To summarise, we take a standard single fluid GRMHD evo¬ 
lution of u'‘,p, Ug, and Pg as a reasonable estimate of the total gas 
properties. This corresponds to assuming that electron conduction 
has a negligible contribution to the dynamics of the total gas and 
that the adiabatic index is independent of the electron thermody¬ 
namic quantities (e.g.; Pglug = [Pp -t Pe\l[ue + «,,] = y - I x 
some function of total gas quantities only). Formally, this assump¬ 
tion requires that the electron internal energy is small compared to 
the proton internal energy. From this, we can calculate the heating 
directly from equation (7) (dropping the conduction terms) without 
requiring an analytic expression for Q. Finally, we use knowledge 
of the nature of heating in a collisionless plasma obtained from 
kinetic theory (described in §5.1) to relate the heating rate per unit 
volume of the electrons, Q,,, to that of the total fluid, Q, and directly 
add it to the electron entropy equation (equation 5), as described in 
§2.2. This completes our model. 

For simplicity, we assume that the adiabatic indices of the 
electron, y^, proton, y^, and total gas, y, are constants, where 
Pk = (jk - l)r<t for k = e,p, or g. This simplifies the numerical im¬ 
plementation of the model, as it allows us to write the entropy per 
particle in a simple form, = (y*-1)“* log(Ptp“”), and avoids the 
complication of having to evaluate TpjTe when updating the total 
fluid variables. This can be seen by noting that 


Pp + Pn 


= (7e - l)(7p - 1) 


1 + TpITp 


i7p - 1) + i7e - DTpITp 


(13) 


which is only constant in the limits that Tp <K Tp or Tp <K Tp. 
From this, we see that this simplification of y = const, is formally 
inconsistent if y^ 7 p, which is generally the case in the accreting 
systems of interest, where the electrons are typically relativistically 
hot (yp X 4/3) but the protons are nonrelativistic (y^ s; 5/3). How¬ 
ever, since equation (13) is bounded between 1/3 and 2/3 and we 
expect Tp < Tp ^ y X yp, we do not anticipate that this approxi¬ 
mation will affect our results significantly. 


2.2 Electron Heating 

We parameterise the heating term, Qp in equation (6) by writing Qp 
= feQ, where fp(J3, Tp, Tp ,....) = QpjQ is the fraction of the total 
dissipation, Q, received by the electrons. This function, in general, 
depends on the local plasma environment and our model is not lim¬ 
ited to any particular choice of fp. As knowledge in the field devel¬ 
ops we can readily incorporate different assumptions about electron 
heating. A more detailed discussion of one physically-motivated 
prescription for fp is given in §5.1. Given a GRMHD solution, the 
total heating rate of a fluid element moving with four-velocity uT in 
the coordinate frame can be computed (from equation 7, dropping 
the conduction terms): 

Q=pTgifdpSg, (14) 

where Sg is the entropy per particle. We can rewrite equation (14) 
in terms of Kg = PgP^^, where Sg = {y— 1)“' log(Kg), as 2 = P^(7 - 


ly^uTdpKg. We use Kg to avoid the undesirable numerical properties 
of logarithms as the argument goes to 0. Likewise, we will often use 
Kp in place of Sp in equation (5). 


2.3 Anisotropic Electron Conduction 


Some care must be taken when generalising the theory of 
anisotropic conduction along magnetic field lines to a relativistic 
and covariant formulation. In particular, the theory must be consis¬ 
tent with causality in that the heat flux should not respond instantly 
to temperature gradients. Our formulation of anisotropic electron 
conduction draws heavily on the treatment of Chandra et al. (2015), 
who consider a single fluid model in which the heat flux is cou¬ 
pled to the dynamics via the stress-energy tensor. We give a brief 
summary of our approach here, highlighting those aspects of our 
electron-only treatment that differ from the formulation in Chandra 
etal. (2015). 

One can derive a perturbation solution for the heat flux, 
by expanding the entropy current in powers of ^ and imposing 
the second law of thermodynamics. The most straightforward rela¬ 
tivistic generalisation of the classical, isotropic heat flux first writ¬ 
ten down by Eckart (1940) is first order in this expansion and 
was later shown by Hiscock & Lindblom (1985) to be uncondi¬ 
tionally unstable, precisely because it violated causality (Chandra 
et al. 2015 showed the same for anisotropic conduction). Israel & 
Stewart (1979) derived a second order solution for which was 
later shown to be conditionally stable (Hiscock & Lindblom 1985; 
Chandra et al. 2015). Here we use a first order reduction of that 
second order model that has been shown to be both stable and 
self-consistent (Andersson & Lopez-Monsalvo 2011). We refer the 
reader to Chandra et al. (2015) for more details. 

We parameterise the heat flux as 

<fp = fif, (15) 

where h'' is a unit vector (hPbp = 1) along the magnetic field four- 
vector, hP, and the scalar f is given by the following evolution equa¬ 
tion: 


((/ipi/) = ( yFs^P’^) = -p 




(16) 


where g is the determinant of the metric, and we used an iden¬ 
tity, VpA^ = dpi -\f^A'^)l to convert covariant derivatives into 
partial ones (eq. 86.9 in Landau & Lifshitz 1975). Here is the 
equilibrium value of the heat flux given by 

0'=“ = -pXe {b^dpTp -t iPapTp). (17) 


where Xe is the thermal diffusion coefficient of the electrons and t 
is the relaxation time scale for the heat flux over which it responds 
to temperature gradients. Note that equation (16) is a relaxation 
equation in which the heat flux relaxes on a timescale t to the equi¬ 
librium value. 

The equilibrium heat flux in equation (17) is the natural rel¬ 
ativistic extension of anisotropic conduction along the magnetic 
field (analogous to the isotropic heat flux of Eckart 1940). The heat 
flux, (fp, then contributes to the electron energy equation as in equa¬ 
tion (5). Physically motivated prescriptions for the parameters Xe 
and T are all that are required to complete the model. We discuss 
one choice of these in §5.1. 


23.1 Stability of Anisotropic Electron Conduction Theory 

In our formalism, we assume that the fluid velocity, m'', and the 
electron number density, Up = plnip, are independent of the elec- 
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tron thermodynamics. Thus, in order to do a perturbative analysis 
we need only perturb the electron temperature, T,,, and the heat 
flux, 0, in equations (5) and (16). Doing this in the fluid rest frame 
in Minkowski space, where m'' = (1,0,0,0), and writing the pertur¬ 
bations in Fourier space as oc exp(/lf ife • T), we find the dispersion 
relation: 


Steps 2 and 3 are described in detail in §3.3 and §3.4, respectively, 
while step 1 will be specific to the choice of the background nu¬ 
merical scheme. 

3.1 Heating in Conservative Codes 


T-+ -+{ye-t)—{b-kf = Q (18) 

r r 

with the solutions: 

4 = ^ |-1 ± y/l - 4(7e - ^)XeT{b ■ . (19) 

The theory is unstable if RefT) > 0, which can only occur if the 
term under the square root is both real and greater than unity. How¬ 
ever, this is impossible for any value of k when > 1, so we con¬ 
clude that equations (5) and (16) are unconditionally stable. This is 
in contrast to the case where equations (5) and (16) are coupled to 
the ideal MHD equations, which is unstable to small perturbations 
unless the relaxation time is larger than a critical value (Hiscock & 
Lindblom 1985; Chandra et al. 2015). 


3 NUMERICAL IMPLEMENTATION OE ELECTRON 
HEATING AND CONDUCTION 


The method outlined above can, in general, be applied to any 
GRMHD “background” simulation. For the rest of this work, how¬ 
ever, we will consider only conservative codes, as the equations of 
ideal MHD can be naturally written in that form. Because of this, 
in what follows we will seek to put all of our evolution equations 
in a conservative form, namely; 


au aF 

dl dx' 


( 20 ) 


where t/ is a “conserved” variable, F' is the corresponding flux in 
the ith direction, and S is the source, which in general includes the 
contribution from the connection coefficients. Equation (20) can 
then be approximated in one spatial dimension by the following 
discretisation: 

(/"+* = U" 




• At 


j+l/2 


•F" 


.0+1/2 


7-1/2 


-. 0 + 1/2 


Ax 


( 21 ) 


where the fluxes are evaluated at face centres using the chosen Rie- 
mann solver. The generalisation to higher dimensions is straight¬ 
forward. 

With that in mind, we can rewrite equation (5): 


[feQ - - ( 22 ) 

where we have used the definition = exp[(yg - 1)^^]. Note 
that equation (22) is a quasi-conservative equation with = 
^E^pu'Ke and - xf^pu'Ke (‘quasi’ conservative because the 
standard definition of conservative equations excludes source terms 
with derivatives). To solve equation (22), we use operator splitting 
in the following sequence of steps: 


1. Solve the conservative equation with =0. 

2. Explicitly update with the heating term (the first term in 
the brackets in eq. 22). 

3. Implicitly solve a matrix equation to include the conduction 
source terms (the rest of the terms in square brackets in eq. 22). 


Formally, the equations of ideal MHD used by conservative 
GRMHD simulations imply that the heating rate per unit vol¬ 
ume, Q, in equation (8) is identically zero. However, conserva¬ 
tive codes implicitly add numerical viscosity and resistivity terms 
to the stress-energy tensor and Maxwell’s equations, respectively. 
The former implies that the numerically evolved stress tensor is in 
fact = Tg'' + + O (truncation error) for some 

numerical viscosity tensor dg', while the latter implies = 0 
-I- O (truncation error). The numerical resistivity can be thought of 
as implicitly introducing a form of Ohm’s law that allows for a 
nonzero electric field four-vector, e'‘. Thus, even though the energy 
implied by + r/M.num Conserved to machine 

precision (see below for details), experiences truncation- 

level heating. This manifests itself as entropy generation: trunca¬ 
tion errors lead to dissipation of magnetic and kinetic energy close 
to the grid scale that is captured as internal energy. We use this 
change in entropy to directly calculate the heating rate per unit vol¬ 
ume of the gas, Q. 

Although is conserved to machine precision, the sec¬ 

ond law of thermodynamics is satisfied only to truncation error. 
Thus there can be locally regions with 2 < 0. In particular, the 
truncation error can be positive or negative, so in places with small 
actual change in entropy or large truncation error the change in 
entropy can be negative. This is the case even in test problems in 
which our methods of calculating the heating give the correct, con¬ 
verged, answer for the fluid variables (see § 4). Thus, while Q may 
be instantaneously or locally negative, when integrated over a suf¬ 
ficient length of time and/or space in the fluid frame it will satisfy 
the second law of thermodynamics. 

We choose this method of calculating the heating rate as op¬ 
posed to introducing an explicit functional form for Q because it 
seems reasonable to assume that for several applications, the grid- 
scale dissipation in conservative codes is a well-defined quantity 
determined by the converged large-scale physics of the problem. 
Turbulence, for example, takes kinetic and/or magnetic energy at 
the largest scales and cascades it down to a small dissipative scale 
where it is converted into internal energy. For a numerical scheme 
with no explicit viscosity, the scale at which dissipation occurs de¬ 
pends entirely on the resolution of the simulation, but we expect 
the heating rate itself will be fixed (in an averaged sense; see, e.g., 
Davis, Stone & Pessah 2010). We expect the same for forced recon¬ 
nection at high ft values, as in the disc midplane, where the large- 
scale dynamics sets the rate at which the field lines of opposite sign 
are brought together. 

The above argument relies on the conservation of energy. In an 
arbitrary space-time, however, conservation of energy is only well- 
defined if the metric is stationary (time-independent) and therefore 
possesses a time-like Killing vector, K>^. If such a vector exists (as 
it does for the Kerr metric of interest in this work), we can construct 
a conserved current from the stress-energy tensor via 7^ = 
where this current satisfies = 0. This allows us to define a 
conserved energy in a coordinate basis: 

E - J' f dx^dx' , (23) 

where the integral is over all space (i.e. the space orthogonal to the 
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time coordinate). Often, the Killing vector takes the form K = d,, 
which simplifies equation (23) to 

£■ = J'-Tl (24) 

Thus, —Tl can be thought of as the conserved energy per unit vol¬ 
ume for a particular choice of coordinates. The total energy, E, is 
conserved to machine precision, modulo fluxes of energy through 
the boundaries, so entropy can only be generated by conversion of 
one form of energy to another. 


3.2 Calculating the Total Heating Rate 

To calculate the heating generated at each time step, we introduce 
an entropy-conserving equation as a reference to compare with 
the energy conservation equation. The entropy-conserving equa¬ 
tion is simply a conservation equation like equation (20) with 
F'k^ = ^Kg = 0- If we call the 

solution to this equation kg (the " denotes the solution correspond¬ 
ing to entropy conservation), then we show in Appendix B1 that 
the total heating rate of the fluid (measured in the fluid rest frame) 
incurred over an interval Ar (measured in the coordinate frame). 


/ p^-‘ 

pu'(Kg - kg) 

Ir-iJ 

At 


where u' accounts for the transformation of At from the coordi¬ 
nate frame to the fluid rest frame, n and n + I denote the values 
at the beginning and the end of the time step, respectively, so that 
4+1 = 4 -I- At, and n -l- 1/2 denotes the values when calculating the 
fluxes. To compute the dissipation rate via eq. (25), we set 
at the beginning of each time step and use when 

calculating the fluxes. Physically, equation (25) means that the La- 
grangian heating rate is set by the difference between the entropy 
implied by the total energy conserving solution (Kg) and the entropy 
implied by the entropy conserving solution (kg). 


3.3 Electron Heating Update 


Let us call k^ the solution to equation (22) without any source terms. 
On top of this adiabatic evolution, electrons receive a fraction, f^, 
of the heating of the gas, = f^Q. In discrete form, this can be 
written as follows. 


(p7ey+l/2 

7e- 1 


(A, -HO""* 


j-n+ljl 


(pA 


,r\n+l/2 


r-1 


-(Kg-kgr 


(26) 


Therefore, the heating update to the electrons, takes 

the following form: 


C' = r ‘ + 21 




(27) 


3.4 Electron Conduction Update 

We note that the evolution equation for the heat flux tj) (equation 16) 
is already in a quasi-conservative form if we define 
and F'^ - We treat the evolution of (j) in an operator split 

way similar to the evolution of Ke with the following series of steps: 

1 Solve the conservative equation with 5^ = 0. 

2 Implicitly solve a matrix equation to include the source terms. 


The source terms in the electron entropy and tj) equation due to 
conduction are given by: 


5*„cond = XpgGe - l)p‘ 

<l> + PXe [b^d^Tg + 


-pyFs 


(28) 


which are discretised in space by using slope-limited derivatives 
across three grid cells and discretised in time by centring the time 
derivatives at 4 + 1 / 2 . The latter discretisation gives us an implicit 
equation for the variables and (p at time f„+i. If we call 0"+' the 
heat flux after being updated by step 1, and ai^ h the electron entropy 
after being updated by heating, then this matrix equation takes the 
form: 


flu 

«21 


ai2 

fl22 


^n+\ 


(29) 


with components 


an = 


/ xFsP^' 

\ Af 


ai2 = [^FgGe- 1)P‘ 


n+1/2 


^21 


xFSP^b'Xe 


«+l/2 




b' 

At 

(1+1 


(30) 


At 


fl22 - 


I xFspit 


\ Af 


and 


bi 


T^gK^xpu' 


At 


■ [(7e - 1)P‘ 


-lrt+1/2 




0+1/2 


/ yFg^ptt' 

\ At 


^ ! - Ip^Xe 




0 + 1/2 


(31) 


0 + 1/2 


0+1/2 


ao+ 1/2 /i.\" . xi 

The system of equations has a straightforward solution. 


^rt+l 


1 


_ biOn - bia2i\ 

011(222 - a21«12 \blCl22 “ I'2«12/ 


(32) 


To ensure that the heat flux, p, does not reach unphysi- 
cally large values, we apply a limiting scheme to keep \p\ < 
(ue -I- pgFjVt.e = <Smax, where pe = pifiglttip and is the electron 
thermal speed. Since we are considering physical systems in which 
the electrons are always at least mildly relativistic, this limit effec¬ 
tively reduces to \p\ < UgCl V3, which corresponds to a ‘saturated’ 
heat flux in which the heat is redistributed at the electron thermal 
speed. The numerical implementation of this limit is to replace the 
values of the thermal diffusivity, Xe, and the relaxation time-scale, 
r, with ‘effective’ values (Chandra et al. 2015): 


Tetf — ++/ 


pm 


(33) 


and 


reff = r/(f^ , (34) 

, <^max / 

© 2015 RAS, MNRAS 000, 1-25 





























Electron Thermodynamics in Black Hole Accretion 1 


where 


m = 1 - 


1 + exp 


1 



+ e, 


(35) 


which sharply transitions from 1 —> e for some small e as |0| 
‘/’max- Thus, according to equation (16), when |(^| > </>max, l</’l decays 
exponentially on a timescale er until it drops below 0max- The 
parameter e is chosen such that the criterion for numerical stability 
is always satisfied (see §3.4.1 and Appendix B3). 


3.4.1 Numerical Stability of Electron Conduction 

A detailed derivation of the criteria for numerical stability is in Ap¬ 
pendix B3. The basic result is that for a Courant-Friedrichs-Lewy 
(CFL) number, C, reasonably chosen between 0 and 1, the relax¬ 
ation time, r, must satisfy 

r>/(C)(^)^„ (36) 

where /(C) is a function of the CFL number. This can be under¬ 
stood as a requirement that the relaxation time t (which we are free 
to choose as arbitrarily large, though which should correspond to a 
physical time scale), must be larger than the time step At (which is 
limited by computational expense) by the ratio between At and the 
standard Courant limit for a diffusive process Afdiff = Ax^lxe- 


3.5 Treatment of the Floors 

Conservative codes deal poorly with vacua of internal energy and 
density. Because of this, many schemes employ floors on internal 
energy and density to ensure that the errors in solving for the prim¬ 
itive variables from the conservative variables do not produce un- 
physically small or negative values. The nature of the model out¬ 
lined above requires special care to be taken when these floors are 
activated, as they introduce artificial changes in internal energy, 
which act as a source of heat, and density, which change the con¬ 
version between entropy and internal energy. 


3.5.1 Electron Energy Floors 

Though the second law of thermodynamics states that the heating 
term u^^d^Kg should be positive definite, numerically we find that 
u'^d^Kg can be locally negative because of truncation errors. This 
introduces the possibility of the electron internal energy going to 
zero (or even becoming negative) due to truncation error fluctua¬ 
tions in our heating term. To correct for this, we implement a floor 
on the electron internal energy that is 1% of the floor on the total 
gas internal energy. That is, if u,, drops below 0.0 In^, we reset u^ to 
0.0 Imj. 


3.5.2 Total Gas Internal Energy Floors 

When the floor on internal energy of the total gas is activated, there 
is an artificial increase in Ug which then shows up in our heating 
term. We treat this addition of energy as if it were a physical, iso- 
choric addition to the energy of the gas and add it to the electrons as 
described above. We emphasise that the internal energy floor does 
not affect the system dynamics in any significant way because it is 
only activated in magnetically-dominated regions where the value 
of the internal energy is dynamically irrelevant. 


3.5.3 Density Floors 

When the floor on density is activated, the total gas internal en¬ 
ergy remains unchanged. Flowever, the value of Ug = kgP^Ky - 1) 
increases by a factor of (fiaooi/PmitY, where pinit is the pre-floor den¬ 
sity. To correct for this, we require conservation of the evolved gas 
entropy when the density floor is activated by decreasing kg by a 
factor of (pfloor/Pinit)^. Furthermore, we enforce that the evolved 
electron entropy remains unchanged by the density floor in the 
same manner by decreasing Ke by a factor of (pfloor/Pimt)’''. Simi¬ 
lar to the internal energy floors, the density floors do not affect the 
dynamics of the system. 


4 TESTS OF NUMERICAL IMPLEMENTATION 

In this section we describe a series of tests that demonstrate the 
robustness and accuracy of our method of evolving the electron in¬ 
ternal energy. We implemented the model described in §2 and §3 
into the conservative GRMHD code, HARM2D (High-Accuracy 
Relativistic Magnetohydrodynamics; Gammie, McKinney & Toth 
2003; Noble et al. 2006). To speed up the computations, we paral¬ 
lelised the code using OpenMP and MPI via domain decomposi¬ 
tion. 


4.1 Tests of Electron Heating 

In what follows we demonstrate the validity and convergence of our 
implementation of electron heating using a number of tests. The 
2nd order convergence of HARM in smooth flows and 1st order 
convergence in discontinuous flows is well documented in Gam¬ 
mie, McKinney & Toth (2003) and we will not reproduce it here. 


4.1.1 Explicit Heating in a Hubble-Type Flow 


To test whether our discretizations of the heating is correctly time 
centred, i.e., converges at the expected 2nd order in time, we fo¬ 
cus here on solving the electron equation when we introduce an 
explicit heating term to the total energy equation. We do this in 
an unmagnetised, ID Hubble-type flow with v oc x (restricting 
ourselves to non-relativistic velocities). In the local rest frame of 
a fluid element, this velocity field gives an outflow in both direc¬ 
tions that is homogenous and isotropic, causing the density to uni¬ 
formly decrease with time as matter leaves the computational do¬ 
main. The velocity profile also scales with time to satisfy the mo¬ 
mentum equation (dv/dt vdvjdx = 0). In the absence of heating, 
the internal energy and pressure evolve according to entropy con¬ 
servation iP oc p^), so that the solution at later times is given by 
(Tchekhovskoy, McKinney & Narayan 2007); 


Vox 

V = - 

1 -I- Vot 

_ tlgfi 

(1+ Votf 
Po 

p = -. 

1 -t- Vot 


(37) 


If we now add a cooling term to the energy equation, of the form: 


_ Ug.oVp (y - 2) 

A+Votf 

the internal energy should evolve as 

^g,Q 

Ug = -T. 

(1 +Votf 


(38) 


(39) 


© 2015 RAS, MNRAS 000, 1-25 









8 S. M. Ressler, A. Tchekhovskoy, E. Quataert, M. Chandra, C. F. Gammie 



Figure 1. LI norm of the error in the electron entropy for heating in a ID 
Hubble-type flow (see §4.1.1). Above a resolution of ~ 1000 the relativis¬ 
tic errors in the analytic result are comparable to the numerical truncation 
errors, so convergence is no longer seen. 


Plugging these solutions in the electron entropy equation (22) for 
fe = l and Ue(t = 0 ) = uq, we obtain: 

«, . ,40) 

7.-2 pJ'(l+«oO 

For the numerical test, we set these analytic solutions as 
the boundary and initial conditions in a one-dimensional grid 
and check if we maintain this solution after a dynamical time of 
L/max[w(f= 0)]. We set 7 = 5/3, 7 ^ = 4/3, max/uoJc) = lO^^c, and 
max(fiVox/Ug) = 1, on a computational domain of 0 < x < 1. For¬ 
mally, since 0^ = kTglmed <K 1, the choice of 7 ^ = 4/3 is unphys¬ 
ical. Flowever the motivation for this choice stems from the fact 
that our primary application is to the inner regions of an accretion 
disc around a hlack hole, where we expect yg x 4/3 t y ~ 5/3. We 
find that our calculation converges at second order (see Figure 1), 
up until the point at which the errors in the analytic solution due to 
relativistic effects become important (which, for m^x(vox) = 10 “^ 
is6Kg/Kg ~ v'^/d ~ 10 ^'’). 

4.1.2 ID (Noh) Shock Test 

In Appendix B4, we show that for a high Mach number shock in 
which the electrons are assumed to receive a constant fraction fg of 
the ‘viscous’ heating in the shock, the post shock electron internal 
energy Ug is given by: 


J 


7+1 

7- 1 


+ 1 + 


7‘ 


1 


— 1 


(41) 


where Ug and 7 are the post-shock internal energy and the adiabatic 
index of the fluid. Equation (41) assumes that the electrons do not 
back react on the shock structure, consistent with the model devel¬ 
oped in this paper. When 7 = y^, equation (41) is equal to //, while 
for 7 = 5/3 and yg = 4/3 it is ~ 0.76fg. Here we check whether 
our numerical implementation of electron heating is consistent with 
this result. 

The initial conditions for this test are an unmagnetised, non- 
relativistic (7 = 5/3), uniform density and internal energy fluid. 
The velocity profile is discontinuous at the center of the grid, with 
a left and a right state given by Vi = -v, = constant > 0. The 
resulting solution is two shocks propagating outwards with a static 
region in between. We focus on a non-relativistic (|n| = lO^^c) flow 
of initially cold gas so that the Mach number of the flow satisfies 


o 

LU 

O 

E 



Resolution, N 

Figure 3. Convergence of the post-shock electron internal energy in the ID 
shock test to the analytic solution (equation 41). The shock’s Mach number 
is ~ 49 (see Sec. 4.1.2). The y. = 7 = 5/3 electrons converge at 1st order, 
as expected, but the yg = 4/3 electrons do not converge to the correct solu¬ 
tion to better than ~ 3% (see Figure 2). This is because our calculation of 
the heating requires a well-resolved shock structure, which is not the case 
for modern shock-capturing conservative codes (see §4.1.2 for details). In¬ 
troducing an explicit bulk viscosity to resolve the shock structure leads to 
convergence for y. 4 7 (see Appendix C). For = 7 , a convenient cancel¬ 
lation makes the evolution of the electron entropy independent of the shock 
structure. 


M » 1. For this test we fix fg = 0.5 and show the results for both 
7 . = 4/3 and yg - 5/3. 

Figure 2 shows the density and electron internal energy as a 
function of position for a shock with M ~ 49 at t = 0.6L/|w;|, 
where L is the size of the computational domain. Figure 3 shows 
that our simulation converges at 1 st order to the analytic result for 
the post-shock electron internal energy when 7 = 7 . but to a value 
differing from the expected result by 3% when yg t y (Fig¬ 
ure 2). This difference is smaller at lower Mach number, as shown 
explicitly in Figure 4. The modest discrepancy between the ana¬ 
lytic post-shock electron temperature and the HARM solution is 
because an accurate, converged calculation of the heating term re¬ 
quires a well-resolved shock structure that gets better resolved at 
higher resolution. This is not the case for modem shock capturing 
techniques, for which the numerical width of the shock is always a 
few grid points and our heating calculation is never able to resolve 
the shock. This is not an issue when 7 = 7 . because the factors of 
density in equation (27) cancel, removing the dependence on the 
shock structure. We show in Appendix C that introducing an ex¬ 
plicit bulk viscosity leads to convergence to the analytical result at 
2nd order for y yg. However, the <3% error as seen in Figures 2 
and 4 is sufficient for our purposes so we do not include a bulk 
viscosity in our calculations. 


4.1.3 2D Forced MHD Turbulence Test 

Another test problem with a known, converged heating rate is 
driven turbulence in a periodic box. If we inject the fluid with a 
constant energy input rate of iJin at large-scales, we should find that 
f QdV - Ein after saturation of kinetic and magnetic energies has 
been reached. Thus, for the electron heating model outlined above, 
after this saturation point the electrons should receive a fraction fg 
of Furthermore, if we have a periodic box in which the par- 
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Figure 2. High Mach number shock results for an electron heating fraction = 0.5 at a resolution of 2000 cells. Top: solid blue line shows the fluid density 
in a numerical simulation. Density undergoes a jump of P 2 /P 1 = (y + Ij/Cy - 1) = 4 at the two shocks, located at x a; 0.35 and x = 0.65. Left: electron 
internal energy relative to total fluid internal energy for = 4/3. The analytic solution is shown with the solid red line and the numerical solution with the 
dotted black line. Right: the same for = 5/3. The analytic solution uses the functional form for Uf/Ug{p) (see Appendices B4 and C for details) and applies 
it to the density returned by the simulation. At this resolution all the fluid variables are essentially converged. The y^ =5/3 electrons show convergence to the 
expected result of Ug = f^Ug (the numerical and analytical lines are essentially on top of each other) while the y^ = 4/3 electrons converge to a value that is 
greater than the analytic result (w^ = 0.379«j for = 0.5; equation 41) by ~ 3%. This is because the internal shock structure is never well resolved without 
an explicit bulk viscosity (see §4.1.2 and Appendix C for details). 


tide number is fixed, then the total internal energy change from 
adiabatic expansion/compression will sum to zero. Thus the an¬ 
alytic result we expect for our model of electron heating is that 
f pTgS^dV = fe f QdV = /c£in. This test checks whether our 
model satisfies this result numerically. 

We start with a static, uniform density fluid with (5 = 6 and 
sound speed Cj,o = 8.6 x lO^^'c in a 2D periodic box. The initial 
magnetic field is uniform: the magnetic field lines are straight and 
lie in the plane of the simulation. Then, at each time step, we give 
Gaussian random kicks to the velocity such that the wave num¬ 
ber satisfies k ■ Sd = 0 (i.e. the driving is incompressible), and 
crl oc k® exp(-8k/kpeak) (compare to Lemaster & Stone 2009). The 
normalisation is fixed such that the rate of energy injection is equal 
to = 0.5pc^g. This leads to a rms turbulent velocity that is 
~ 0.8cs,o ~ 1.8 ua, so that the turbulence is subsonic and roughly 


Alfvenic. The peak driving wave number is set to half the box size: 
fcpeak = 47r/L. Furthermore, we ensure that no net momentum is 
added to the box by subtracting off (from the kicks) any average 
velocity that would have been generated by the kicks. For this test 
we fix fe = 0.5, y = 5/3 and y^ = 4/3. 

Figure 5 shows our results for the electron and total internal 
energies as a function of time in the box at 512^. We see that once 
approximate saturation of the turbulence is reached at f ~ L/Ci,o (or 
f ~ 0 in the figure), the internal energies are in very good agree¬ 
ment with a linear fit, as expected given the constant rate of en¬ 
ergy injection. For a parameterization of f u^dV = get + be and 
J UgdV = ggt + bg, ge and gg represent the electron and total heating 
rates, respectively. These can be compared to the energy injection 
rate, which is a fixed constant of O.Spd^g. At a resolution of 
512^ , we find the total heating rate differs from by ~ 4%, while 
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Figure 4. Percent error in the post-shock electron internal energy for jg = 
4/3 and y = 5/3 as a function of Mach number in the ID shock test as 
computed by HARM at a resolution of 2000 for a fractional heat given to 
the electrons of/(, = 1/2 and an initial M(,/u^ = 0.1. The analytic solution is 
given by equation (B27). The final time was fixed such that the two shocks 
were located at x = 0.25 and x = 0.75 in a 0 < x < 1 domain. Note that 
the fractional errors are always g 3.3%. The change in the percent error as 
Mach number goes to 1 is because the flow becomes increasingly smooth 
and the electron internal energy is no longer converged at a resolution of 
2000 . 


Table 1. Turbulence Test Linear Fits (§4.1.3) 


Resolution: 

128 

256 

512 

gg - Q.SEiA 

0.0027£',„ 

0.0042F,-,, 

0.0024£:,„ 

Sg ~ 

-0.000 

0 .0012F,„ 

-0.0016£:,>, 


Fractional error in the electron heating rate relative 
to the analytic solution. 

^ Fractional error in the total heating rate relative to 
the analytic solution. 


the electron heating rate differs from fgEi„ hy ~ 2%. Unfortunately, 
a rigorous convergence study of these quantities is not possible be¬ 
cause of the nature of turbulence in 2D. Due to the inverse energy 
cascade, the kinetic and magnetic energies never truly saturate and 
convergence of any of the fluid variables is never achieved. This 
can be seen from Table 1, where the values of gg and gg are quoted 
at various resolutions, neither of which display significant conver¬ 
gence to the expected 0.5£,„ and £,,,, respectively. Nevertheless, 
we find the percent level error found at all resolutions to be suffi¬ 
ciently small to satisfy our error tolerance in the full accretion disc 
simulation. 


4.1.4 Shadow Solution 

For our two-temperature model, we can seek a solution in which 
the electron fluid simply ‘shadows’ the total gas, in that Ug oc Ug. 
Such a solution is found by setting Ug{t = 0) = fgUg(t = 0) at some 
initial time, since we can solve the first law of thermodynamics for 
the electrons with Ug = fgUg for all time, assuming that jg = y and 



Figure 5. Electron and total gas internal energies summed over the grid as a 
function of time at a resolution of 512^ for the forced subsonic MHD turbu¬ 
lence simulation with an assumed electron heating fraction of fg = 0.5 and 
initial /J = 10. We define t = {t - ti)CgfllL andUg fg) = Ug^^g) — « 5 ,(c)(t = t;), 
where f, ~ ALjCsf) is the time at which kinetic and magnetic energy roughly 
saturate. We normalised the integrated internal energy by the energy injec¬ 
tion rate, = 0.5pd q , so that the y-axis has dimensions of time which 
we measure in units of Llcgp. In these variables, the analytic solutions for 
the total gas and electron internal energies are lines with slopes of I and 
fg = 0.5, which are plotted as solid lines, to be compared to the simulation 
results which are represented by points. We find that the electron heating 
rate is 0.5 of the total heating rate, consistent with the analytic solution 
given the input value of fg = 0.5. For a numerical comparison at differ¬ 
ent resolutions, see Table 1, which shows the results of applying a linear 
regression fit to the internal energies. 


fg is a constant. This can be seen from the first law: 

i/df^Ug = fgpTgi/df^Sg - - - -u'^dfip, (42) 

because the electrons will always get a fraction, fg of the entropy¬ 
generated heat (the first term on the RHS of equation 42), while the 
compression term is directly proportional to Ug oc fgUg. This solu¬ 
tion is valid regardless of the details of the overall fluid evolution, 
so we can apply it to an arbitrarily complicated system. 

For this test, we evolve the electron internal energy in the full 
accretion disc simulation around a rotating black hole as outlined 
in §5. We initially apply small {Sugjug ~ 0.04) perturbations to the 
electron internal energy about the average value of Ug^gp = O.Sugp, 
and set fg = 0.5 and jg = y = 513. For this test alone, we set 
the floor on electron internal energy to be a fraction fg of the floor 
on the total fluid (as opposed to our usual choice of 1%). If this 
latter step were neglected, the floors would cause the polar regions 
to differ significantly from the expected result, though leaving the 
disc and corona unaffected (i.e. they still satisfy the analytic result). 
The test is whether or not our simulation can maintain this result 
over the run time of 2000M. 

Running this test at a resolution of 512^ gives an average frac¬ 
tional error of 


, N-l JV-1 
;=0 1=0 


([“c/Ms].y “ /<■) 


0.8%, 


(43) 


which is smaller than our initial perturbations and shows that our 
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numerical solution correctly evolves equation (42) even in a com¬ 
plex problem with MHD turbulence, weak shocks, and other heat¬ 
ing processes in the presence of a curved metric. 

4.2 Tests of Electron Conduction 

Our model and testing suite for conduction closely resembles that 
of Chandra, Foucart & Gammie (2015), so we leave the details to 
Appendix A. In summary, we have found second order convergence 
for linear modes, for a static, ID atmosphere in the Schwarzschild 
metric, and for a relativistic, spherically symmetric Bondi accretion 
flow. We also show that the electrons properly conduct along field 
lines in a 2D test. 


5 APPLICATION TO AN ACCRETING BLACK HOLE IN 
2D GRMHD SIMULATIONS 

We apply the new methods discussed in §2 and the numerical im¬ 
plementation described in §3 to the astrophysical environment of 
an accretion disc surrounding a spinning black hole as described 
by the Kerr Metric with a spin parameter of a = 0.9375. For this 
spin the last stable circular orbit is sj 2.04rj and the thin disc ra¬ 
diative efficiency is =» 0.18 (Novikov & Thome 1973). We use the 
conservative code HARM (Gammie, McKinney & Toth 2003) as 
our background GRMHD solution. Our initial conditions for the 
total fluid are the Fishbone & Moncrief (1976) equilibrium torus 
solution (see Appendix D) with inner radius Hn = 6rj and with 
the maximum density of the disc occurring at r^ax = I2rg. Note 
that here and throughout r and 9 refer to the Boyer-Lindquist coor¬ 
dinates. This equilibrium solution has a temperature maximum of 
7.5 X lO'^K and a thickness"* of hjr ~ 0.18 at r^ax- We normalise 
the torus density distribution such that the maximum value of den¬ 
sity in the torus is Pmaxd = 1 and perturb the internal energy of the 
gas with random kicks on the order of 6 ug/ug ~ 0.04 to provide the 
perturbations for the magnetorotational instability (MRI, Balbus & 
Hawley 1991) to develop.^ We overlay this equilibrium solution 
with an initial magnetic field with 2f’n,ax/^max = 100 (where max 
refers to the maximum value inside the torus), defined by the scalar 
vector potential: 

oc (p/Pmax - 0.2) cos 9, (44) 

if p > 0.2pn,ax and 0 otherwise. This vector potential defines 
two meridional loops contained in the toms that are antisymmet¬ 
ric about the equator. This choice ensures that the field lines are 
not along constant density. Since constant density implies constant 
temperature when entropy is constant, field lines along constant 
density would be isothermal in the initial condition (as would hap¬ 
pen if we dropped the factor of cos0 in eq. 44). 2D MHD torus 
simulations are unable to reach a statistical steady state in which 
the initial conditions are forgotten, so initially isothermal field lines 
could artificially suppress electron conduction even at later times. 
We choose the 2-loop initial condition to avoid this. 

"* Here we define fi/r = Jf pu'\d — 7i/2\ x/^ddd^j jjpu' . 

^ In addition to the electron specific floor described in §3.5.1, there aie 
also floors on the density and internal energy of the HARM single fluid 
GRMHD solution. These are pnoorC^ = max [fe^/50,10^‘*(r/rg)“^''^pn,axC^] 
and Hfioor = max [&^/250,10“^(r/rj)“ "^/^PniaxC^]. Note that the unit choice 
for the background ’’atmosphere” is such that the initial torus maximum 
density is PmaxC^ = 1 and the initial torus internal energy is Umax ~ 0.01. 


For the electrons, we start with Uj/Ug = 0.1, and run two dif¬ 
ferent models for /<,, described below. For conduction runs, we set 
the initial heat flux to zero. 

Our conduction runs are all in 256“ grids with a physical 
size of the domain in spherical polar coordinates of (Ri„,Rout) x 
(^?in>^^out) = (0.8rH, lOOOrj) x (0,;r), where rjj = rg(l + Vl - a^) is 
the black hole event horizon radius. Fora = 0.9375, th ~ 1.35r^. In 
the regions with r < 50r^, the code uses modified Kerr-Schild co¬ 
ordinates (t, v', x^, and ip) of Gammie, McKinney & Toth (2003), 
so that the regions with the highest resolution are near the mid¬ 
plane close to the horizon. For r > 50r^, we use hyper-exponential 
coordinates to move out the outer radial boundary r = and 
limit unphysical reflection effects by defining the internal code co¬ 
ordinate x* implicitly by the equation (Tchekhovskoy, Narayan & 
McKinney 2011): 

exp(v*) : r < 50rg 

expjx*-I-[x* - x'(r = fiOrg)]"*} : r > 50rg. 

The electron heating-only (i.e. without conduction) runs have the 
same parameters but a higher resolution of 512^. At the inner 
and outer radial boundaries we apply the standard outflow (copy) 
boundary conditions, at the polar boundaries we apply the standard 
antisymmetric boundary conditions (with all quantities symmetric 
across the polar axis except u® and B®, whose signs are reversed). 

Figure 6 shows the background HARM solution for the den¬ 
sity, magnetic field, temperature, plasma p = 2 PglE, and the heat¬ 
ing rate per unit volume in the coordinate frame, -Qu,, averaged 
over the time interval 900 - 1100 Xg/c, as well as the initial field 
configuration. After ~ 1200rg/c the turbulence starts to decay, an 
artefact of 2D simulations in which MRI turbulence is not sustain¬ 
able. 

As noted in §3.1, we find locally that Q < 0 (violating the 
second law of thermodynamics) in many regions due to truncation 
errors. This is because HARM satisfies the total energy equation 
to machine precision but only satisfies the second law of thermo¬ 
dynamics to truncation error. However, while Q may be instanta¬ 
neously or locally negative, when integrated over a sufficient length 
of time and/or space in the fluid frame it will satisfy the second law 
of thermodynamics. In our torus simulation, for instance. Figure 7 
shows that when averaged over 9 and time (900 - 1100 Xg/c), the 
heating rate is entirely positive definite within the region of interest. 
Furthermore, when integrated over the volume enclosed between 
the event horizon, rn, and r = 6rg (roughly the radius at which the 
accretion time ~ 1000 Xg/c), we find 

2n K 

-u,Q{d -I- a^ cos^ 6 ) sin 6 dr d0 d(p ~ 0.17mA, (45) 

0 0 r// 

where the factor of -u, converts Q to the coordinate frame. In equa¬ 
tion (45), M is the accretion rate of the black hole in terms of coor¬ 
dinate time (corresponding to time measured by a distant observer) 
at the event horizon radius, r = rn, 

M= J' pu'^(d + cos^ 9) sin 9 d9d<p. (46) 

T=rin 

The heating rate in equation (45) is in excellent agreement with 
that expected for a rapidly spinning black hole (e.g., the Novikov 
& Thome 1973 model predicts a radiative efficiency of » 0.18 for 
a = 0.9375). 

In this work, all mass-weighted averages are computed using 
the weighting function: pu' which represents the conserved 
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mass per unit coordinate volume. For example, a radial average of 
a function f(x\x^, x^,t) is computed as: 


•^max 

f f{x^ ,t)pu’ 

■^min 

•^max 

J pu’ xE^dx^ 

JC* . 
nun 


(47) 


5.1 Electron Parameter Choices 


Here we describe physically motivated estimates of the electron 
heating fraction, fg, and the electron thermal diffusivity, appro¬ 
priate for low-collisionality accretion flows such as that of Sagittar¬ 
ius A*. A more comprehensive exploration of physical models will 
be explored in future work. 

We consider two simple models for the electron heating frac¬ 
tion E. The first sets = 1/8, a constant. Because the electron 
adiabatic index is not the same as the proton (total) adiabatic index, 
and the heating is not spatially uniform, a constant /<, model does 
not necessarily lead to a constant TpIT^. The second, more physi¬ 
cal model, sets fe based on theoretical models of the dissipation of 
MHD turbulence in low-collisionality plasmas. These genetically 
predict that electrons receive most of the turbulent heating at low p 
while protons receive most of the turbulent heating at high p. This 
is true both for reconnection (Numata & Loureiro 2015) and col¬ 
lisionless damping of turbulent fluctuations (Quataert & Gruzinov 
1999). This dependence on p is the key qualitative feature of our 
chosen model of //. For concreteness, we use the specific calcula¬ 
tions of Howes (2010) who provided a simple fitting function for 
the electron to proton heating rate as a function of plasma parame¬ 
ters in calculations of the collisionless damping of turbulent fluctu¬ 
ations in weakly compressible MHD turbulence like that expected 
in accretion discs. These models do a reasonable job of explaining 
the measured proton and electron heating rates in the near-Earth so¬ 
lar wind (Howes 2011). The functional form of is derived from 
the relations: 


n ^2 fl2-0.21ogio(7'p/7'.) 

\cp _ ^2 “T Pp 

a “''‘c2+^2-o.2iog,„(r,/r.) 



tlUp 


(48) 


with Cl = 0.92, C 2 = Lb/fTp/TQ, and Cg = 18 -H 5 logiofTp/TQ for 
TpjT, > 1, while Cg = l.lKTpIT,) and Cj = 18 for T„/T, < 1. The 
corresponding result for fg is simply 


fe = t;- 77 = 1 — TTTTr- 

Qp Qe 1 Qp! Qe 

The critical assumption used in deriving equation (48) is that the 
turbulent fluctuations on the scale of the proton Larmor radius 
have frequencies much lower than the proton cyclotron frequency. 
This is believed to be well-satisfied for weakly compressible MHD 
turbulence in accretion disks (e.g., Quataert 1998). For concrete¬ 
ness, we note that for TpjTg - 1 and Pp - (0.1,0.3,1,10), we 
have QpIQg = (0,0.01,0.16,8.6), while for Tp/T, = 10 and 
Pp = (0.1,0.3,1,10), QpIQg = (0,0.001,0.09,12), respectively. 
This demonstrates the strong transition from predominantly elec¬ 
tron to predominantly proton heating with increasing Pp, with the 
transition happening at a value of Pp that depends weakly on the 
proton to electron temperature ratio. This implies that we expect 
strong electron heating in the corona and jet regions but suppressed 
electron heating in the bulk of the disc. 

We reiterate that the key feature of equation (48) is not the 


precise value of the predicted QpIQe, but rather the transition from 
Qp <: Qe for j8p » 1 to (2p <K Qe for « 1. This qualitative 
transition is much more robust than the specific functional form in 
equation (48) (e.g., Quataert & Gruzinov 1999; Numata & Loureiro 
2015). 

For the electron thermal diffusion parameters, since Xe is a 
diffusion coefficient, we assume that it has the form 


Xe = Oecr, 


(50) 


where Ug is a dimensionless thermal diffusivity, and r is the radial 
distance from the center of the black hole, which is comparable to 
the density scale height of the disc, H. Since we are interested in 
fairly relativistic electrons, we choose the relevant velocity to be 
c in our diffusivity estimate. In what follows, we consider a range 
of dimensionless diffusivities, ~ 0.1 - 10. A typical value of 
ttf ~ 1 is motivated by the idea that particles scatter roughly after 
moving a distance comparable to the length-scale over which the 
magnetic field strength, density, etc. change. In fact, for high beta 
plasmas, the mean free path due to wave-particle scattering can be 
significantly lower, reducing the thermal diffusivity significantly. 
In Appendix B2 we discuss the specific limits imposed by electron 
temperature anisotropy instabilities present in a turbulent plasma. 
In particular, the whistler and firehose instabilities lead to limits on 
IxTelTe (eq. B9) and thus the electron viscosity and thermal diffu- 
sivity, where the temperature anisotropy is defined with respect to 
the local magnetic field. In terms of the electron thermal diffusivity, 
this becomes= min(aerc,Xm!ix), where;^^n,ax is set by velocity 
space instabilities and is estimated in Appendix B2. Finally, we 
choose the relaxation time scale, t to be given by the thermal time 
scale: 


T 




(51) 


Comparing this to the stability condition given by equation (36), 
we see that stability is ensured if 


v,i, < 


Ax 
At' 


(52) 


Since we have a non-uniform grid, the time step Ar is essentially set 
by the light crossing time of the smallest grid cell (i.e. that nearest 
the horizon), meaning that At < cAx near the horizon and At « 
cAx further from the horizon. For a reasonable choice of a CFL 
number of 0.5, we find that equation (52) is satisfied everywhere 
and is not a limiting factor in our simulation. Moreover, we find that 
the exact value of r is not critical as long as it satisfies numerical 
stability and is not too long (e.g. is less than a local dynamical 
time). 


5.2 Electron Heating Only 

In this section we focus solely on the effects of separately evolv¬ 
ing the electron internal energy equation without conduction in our 
black hole torus simulation and compare the results for different 
electron heating models. 


5.2.1 Constant Electron Heating Fraction 

Figure 8 shows the temperature ratio, TgjTg, averaged over the in¬ 
terval 900 - 1 lOOrj/c for fe - I /8. We reiterate that Tg here is the 
temperature inferred from the underlying single fluid GRMHD so¬ 
lution (approximately the proton temperature in our model) while 
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Figure 6. Properties of our 2D black hole accretion simulations. The top panel shows the density over-plotted with magnetic field lines in the initial conditions 
(left) and averaged over 900 - 1100 r^/c (right). The remaining panels are the total gas temperature in units of nipd (middle left), the plasma parameter, 
p = IPgjb^ (middle right), and the absolute value of the heating rate per unit volume in the coordinate frame, \Qu,\, in units of Mc^/{ d~S) (bottom), all 
averaged over time in the interval 900 - 1100 fj/c. Note that for calculating the average yS, we use 2{Pg)/{b^), where () denotes an average over time. These 
plots represent the background GRMHD solution on top of which we separately solve the electron entropy equation. 


Tg is the electron temperature determined from our separate elec¬ 
tron entropy equation. We include this constant result primarily 
because it is conceptually similar (although quantitatively differ¬ 
ent) to the constant Tp/Tg assumption often used in the literature. 
Notice that the resulting TejTg ratio, seen in Fig. 8, is non-uniform 


despite the constant f^. Also note that due to the fact that MHD 
turbulence is unsustainable in 2D simulations, the heating dies off 
affer ~ 1200rj/c of evolution and prevents the outer r > lOr^, re¬ 
gion of the disc from ever being heated substantially. However, as 
we will see later, conduction can occur at a much faster (electron 
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Figure 7. Mass-weighted average (see eq. 47) of the heating rate per unit 
volume in the coordinate frame, averaged over time in the interval 900 - 
1100 Tj/c and over Q from 0 to n. Note that for our metric sign convention, 
u, < 0. The total volume integrating heating out to ~ is ~ 0.17Mc^ 
(equation 45), comparable to the Novikov & Thorne (1973) heating rate for 
this black hole spin. 



Figure 8. Ratio of (Jg) l(Tg) in our black hole accretion simulation, where 
() denotes an average over time in the interval 900 — 1100 r^/c (where 
Tg is the temperature of the single fluid GRMHD simulations and Tg is the 
electron temperature). These results assume a constant fraction of dissipated 
heat is given to the electrons (fg = 1/8). Compare with the more physical 
/1-dependent heating results in Figure 9. 


thermal) speed along the magnetic field lines and can affect the so¬ 
lution at somewhat larger radii. 


5.2.2 P-Dependent Electron Heating 

Figure 9 shows the temperature ratio, TgjTg, the electron tempera¬ 
ture itself, &g, and the electron heating fraction, fg, averaged over 
the interval 900 - 1100 r^/c for the //-dependent heating model of 
equation (48), which we regard as a more physical electron heating 
model than fg = const. We note that this leads to hot electrons being 
strongly concentrated in the corona of the torus in between the disc 
and the jet, where // is the smallest (and fg~l from equation 48). 


This is also clear from the ID profiles of electron temperature as a 
function of polar angle in Figure 10. 

Figure 10 shows the mass-weighted average over radius (r = 
5—lrg) of the electron and gas temperatures, plotted versus the po¬ 
lar angle, 9. The fg = 1/8 electrons and total gas temperatures have 
mild variation in T with 6 , while the fg = fg(fl) electrons have 
significantly higher temperatures in the polar regions. This demon¬ 
strates that the non-uniformity of the electron temperature in the 
feiP) model is primarily caused by the strong //-dependence of our 
model of fg as opposed to any non-uniformity of the heating rate 
itself (Figure 6). 


5.3 Conduction and Electron Heating 


We now consider the effects of electron conduction on the electron 
temperature structure of black hole accretion discs. We focus on 
the more physical model of //-dependent heating described in §5.1. 
In all of our calculations, we include the velocity space instability 
limit on the electron thermal conductivity (Appendix B2), although 
runs without this limit produce similar results because // is mod¬ 
est (< 1 - 10) in the inner regions of these simulations (Figure 6). 
Figure 11 shows the electron temperature as a function of radius at 
the mid-plane in the simulations with and without conduction. Fig¬ 
ure 12 shows the effects of conduction more quantitatively via the 
fractional change in temperature between the electron temperature 
solution with conduction and that without. 

To summarise Figures 11 and 12, conduction has little effect 
on the electron temperature for Ug < 1. However, for Og > I, con¬ 
duction leads to a significant radial redistribution of heat such that 
the electron temperature is factors of a few larger at large radii. 
Even for Ug > I, however, the angular redistribution of heat is much 
less efficient, as seen in the radially and time-averaged electron 
temperatures in Figure 10 for Ug = 10. This is primarily because 
of the structure of the magnetic field, as can be seen by noting that 
the regions where conduction modifies the temperature in Figure 12 
largely follow magnetic field contours which do not efficiently con¬ 
nect the polar and equatorial regions. To aid the interpretation of 
these results. Figure 13 shows the heat flux f normalised to the 
maximum value = dig + pgC^)v,/, even for Ug = 10 the heat 
flux is still well below the saturated value in significant parts of 
the domain. We now summarise and interpret these results in more 
detail. 

For Ofj < 1 we find conduction to have only a small effect 
on the electron thermodynamics in the accretion disc, despite the 
relatively high conductivity. We can understand this result as being 
due to the suppression of the isotropic heat flux by being projected 
along field lines, quantified by the ratio, 




(53) 


where and are evaluated using the electron temperature 
as evolved without conduction and which we now define. For this 
diagnostic, we use 

ql„ = -pXeff''{d,Tg + Tga,), (54) 

where /z''’' = 4- g'^’' is the projection tensor that projects along a 

space-like direction perpendicular to the fluid velocity u''. This pro¬ 
jection ensures that the heat flux in the fluid frame has a zero time- 
component. Likewise, for we use the first order anisotropic 
heat flux: b''. Note that in equation (53), e is always 
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Figure 9. Ratio of (T,,) /{Tg) (top), electron temperature, (T,,), in units of 
irigA (middle), and electron heating fraction, (f,,) (bottom), where () de¬ 
notes an average over time in the interval 900 - 1100 fj/c. These results are 
for yS-dependent heating (see §5.2.2). Compare to Figure 8 for a constant 
electron to proton heating ratio. The highly non-uniform distribution of /? 
(see Figure 6) and the strong /} dependence of the electron-to-total heating 
ratio (equation 48) lead to a strong angular dependence of Tg/Tg. 



Figure 10. Mass-weighted average of total gas and electron temperature 
(in units of nigd) as a function of the polar angle, 6. We show the electron 
temperature with and without conduction for a /3-dependent electron heating 
fraction, fg, as well as without conduction for a constant electron heating 
fraction fg = 1/8. The results are averaged over time from 900 - 1100 
rg/c and averaged over r from 5 -1 Tg. Note that the total gas temperature 
has been multiplied by a constant fraction to more clearly compare to the 
electron temperatures. The electron temperature with /3-dependent heating 
displays much stronger 6 variation because the electron heating fraction 
itself values with 8 (see Figure 9). Conduction has only a modest elfect on 
redistributing heat in 8 due to the geometry of the field. 


< 1 because both heat fluxes are mutually orthogonal to Fig¬ 
ure 14 shows l^anisol / l?isol in ouf torus simulation, where we find 
the suppression of the isotropic heat flux to be around e ~ 0.2. 
The simplest explanation for this small number is that the field is 
predominantly in the (p direction, where the temperature gradient 
is identically 0 in 2D simulations. For instance, in local shearing 
box calculations, Guan et al. (2009) found that the typical angle 
between B and tp was ~ 10 - 15°, corresponding to a suppression 
of the heat flux with e ~ 0.25. 

Contrary to the Og < I cases, setting ag > I causes conduc¬ 
tion to have a significant effect by redistributing the electron heat 
from the coronal regions to the bulk of the torus at larger radii. This 
redistribution of heat causes the electron temperature to actually 
exceed the total gas temperature in certain regions, which formally 
violates our assumption that Tg <K T,,. 

While the calculation with Og = 10, or with Xe = lOrc, 
might seem to use an unphysically large conductivity, roughly cor¬ 
responding to a length scale for conduction of ~ 10//, where H is 
the disc density scale height, the heat flux in these calculations is 
limited to be smaller than the value set by the physically motivated 
whistler criterion in equation (BIO) and to be less than the saturated 
heat flux ~ UgC. As Figure 13 shows, the heat flux is saturated in 
only part of the domain. Furthermore, the appropriate length scale 
for conduction should be the scale height along field lines, which 
could be significantly greater than the overall density scale height 
if the field has a large toroidal component. For these reasons, we 
believe that the larger Ug solutions may in fact be physical because 
they correspond to a heat flux closer to the saturated value ~ UgC 
expected in low-collisionality plasmas. 
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Figure 11. Electron temperature in the mid-plane (8 = n/2) in units of 
for black hole accretion simulations with /3-dependent heating and for elec¬ 
tron conduction with dimensionless conductivity = 0,0.1,1,10 (where 
the electron thermal diffusivity is;^'^ = a^rc; see §5.1). The results are time 
averaged over the interval 900 - llOOr^/c. For alphas ^ 1, conduction re¬ 
distributes energy from small to large radii, increasing the electron temper¬ 
ature at larger radii. Compare to Figure 10, which shows that redistribution 
of heat in the polar direction is less efficient. 


6 CONCLUSIONS 

We have presented a method for evolving a separate electron en¬ 
tropy equation in parallel to the standard equations of ideal General 
Relativistic MHD. Our motivation is the study of two-temperature 
radiatively inefficient accretion flows (RIAFs) onto black holes, in 
which the electron-proton Coulomb collision time is sufficiently 
long that the proton and electron thermodynamics decouple (e.g., 
Rees et al. 1982). Understanding the electron temperature distribu¬ 
tion close to the black hole is necessary for robustly predicting the 
radiation from the numerical simulations of black hole accretion 
(and outflows) in the sub-Eddington regime. 

The long-term goal of the present work is to incorporate the 
key processes that influence the electron thermodynamics in RIAFs 
into GRMHD simulations: heating, thermal conduction, radiative 
cooling, and electron-proton Coulomb collisions. In the present pa¬ 
per we have focused on the first two of these processes. Specifically, 
we have developed, implemented, and tested a model that quanti¬ 
fies the rate of heating in a conservative GRMHD simulation (§2). 
We then assign a fraction of this heating to the electrons based 
on a microphysical model of the key heating processes (e.g., turbu¬ 
lence, reconnection, shocks; see, e.g., §5.1). In addition, we have 
implemented and tested a model of relativistic anisotropic conduc¬ 
tion of heat (by electrons) along magnetic field lines, based on the 
Chandra et al. (2015) formulation of anisotropic relativistic con¬ 
duction (§2.3). The electron thermal diffusivity is a free parameter 
in this calculation. For the black hole accretion disc applications 
of interest, we advocate a ‘saturated’ heat flux in which the ther¬ 
mal diffusivity is rc, subject to additional constraints imposed by 
velocity space instabilities and scattering by wave-particle interac¬ 
tions (Appendix B2). 

We implemented our electron energy model in a conserva¬ 
tive GRMHD code HARM2D (Gammie, McKinney & Toth 2003), 
though the model we have developed can be applied to any un¬ 
derlying GRMHD scheme. For simplicity, the implementation in 
this paper neglects the back reaction of the electron pressure on 



Figure 12. Fractional difference in electron temperature between solutions 
with and without electron conduction shown in colour (see colour bar for 
details) over-plotted with magnetic field lines shown as solid black lines. 
The fractional difference is calculated as / (T^.o)- 1, where () denotes 
an average over time from 900- 1100 r^/c. The results include //-dependent 
electron heating for = 0.1 (top), = I (middle), and = 10 (bottom 
panel), where the electron thermal diffusivity (§5.1). Higher 

allows more heat to flow from the inner regions to larger radii. For =0.1 
conduction has a negligible effect on the electron temperature, while for 
> 1 conduction leads to order unity changes in . 
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Figure 13. <|</>|)/(^^'max )5 the ratio of the electron heat flux to the maximum 
value ^max = (Wg + peC^)V{^e, where () denotes an average over time from 
900 - 1100 Tgjc. This is calculated based on the results of a black hole 
accretion simulation with ^S-dependent electron heating and a dimension¬ 
less electron thermal conductivity of ag = 10 (where the electron thermal 
ditfusivity is;^'^ = see §5.1). Comparison with Figure 12 shows that 
conduction has a significant effect on redistributing heat only in the regions 
where the heat flux is saturated or nearly saturated. However, even for a high 
electron thermal conductivity of ag = 10, the heat flux is still well below 
the saturated value over much of the domain. 
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Figure 14. {l5anisol)/(l‘?isol), the ratio of the anisotropic (field-aligned) heat 
flux to the isotropic heat flux , where () denotes an average over time from 
900-1100 fj/c. This is calculated based on the results of a black hole accre¬ 
tion simulation with yS-dependent electron heating but without conduction. 
The factor of ~ 5 - 10 suppression of the field aligned heat flux is roughly 
consistent with that expected from local shearing box calculations of MRI 
turbulence, where B is aligned with the 0 direction (e.g. Guan et al. 2009). 

the dynamics of the accretion flow. We believe that this is a rea¬ 
sonable first approximation given some of the uncertainties in the 
electron physics. Formally, this approximation is is valid only when 
Tg <K Tp though we expect it to be a reasonable first approximation 
when Tg < Tp in regions with plasma /3 > 1, i.e., in the regions 
where gas thermal pressure forces are dynamically important. 

We have demonstrated that our implementations of electron 
heating and conduction are accurate and second order convergent 
in several smooth test problems (§4 and Appendix A). For shocks. 


the heating converges at first order but to a post shock temperature 
that differs from the analytic solution by < 3% when the electron 
adiabatic index differs from the adiabatic index of the fluid in the 
GRMHD solution (e.g.. Figure 2). This discrepancy arises because 
standard Riemann solvers ‘resolve’ the shock structure with only a 
few grid points. Including an explicit bulk viscosity to broaden and 
resolve the shock leads to a converged numerical solution for the 
post-shock electron energy that agrees with the analytic solution 
(Appendix C). In practice, the < 3% discrepancy between the nu¬ 
merical and analytic solutions for standard Riemann solvers is suf¬ 
ficiently accurate given other uncertainties in the electron physics. 
For this reason, we do not use bulk viscosity in our calculations. 
Moreover, strong shocks are rare and account for a negligible frac¬ 
tion of the dissipation in accretion disc simulations with aligned 
black hole and accretion disc angular momentum. 

In addition to formulating and testing our electron energy 
equation model, we have also presented a preliminary applica¬ 
tion of these new methods to simulations of black hole accretion. 
Specifically, we have studied the impact of realistic electron heat¬ 
ing and electron thermal conduction on the spatial distribution of 
the electron temperature in 2D (axisymmetric) simulations of black 
hole accretion onto a rotating black hole. We find that the resulting 
electron temperatures differ significantly from the assumption of a 
constant electron to proton temperature ratio used in previous work 
to predict the emission from GRMHD simulations (Moscibrodzka 
et al. 2009; Dibi et al. 2012; Drappeau et al. 2013); see, e.g. Fig¬ 
ures 9-11. This is due to the strong /(-dependence of the electron 
heating fraction, fg, described in §5.1: electrons are preferentially 
heated in regions of lower ft, causing TgfTp to be larger in the coro¬ 
nal regions compared to the midplane. In addition, we find thaf fhe 
effect of thermal conduction on the electron temperatures is sup¬ 
pressed by the fact that the heat flux must travel along field lines, 
which are predominantly toroidal and thus not aligned with the 
temperature gradient. Specifically, we find that electron conduction 
modifies the temperature distribution only if the effective electron 
mean free path along the magnetic field is > the local radius in the 
flow (see Figure 12). In this case, there is a net transfer of heat 
from the corona to the bulk of the disc. This increases the electron 
temperature at larger radii by a factor of ^ 2. 

It is important to stress that the unsustainability of MHD tur¬ 
bulence in 2D simulations (e.g., Guan & Gammie 2008) limits how 
thoroughly we can interpret the accretion disc results presented in 
this paper. Since a steady state is never truly reached, the bulk of 
the disc retains memory of the initial conditions and only the inner¬ 
most regions (r < lOr^) develop significant turbulence. This could 
artificially limit the effects of electron conduction because the ther¬ 
mal time for relativistic electrons is ~ rjc and is thus substantially 
shorter than the local dynamical time only at large radii. Future 
work will use the methods developed here in 3D simulations. 

It is also important to stress that, as in previous work, our re¬ 
sults for both the gas and electron temperature are not reliable when 
E » pd. In these regions the ratios of Efpd and EfUg are so 
large that the evolution of the density and internal energy are dom¬ 
inated by truncation errors in the magnetic field, to which they are 
nonlinearly coupled by the total energy equation. This requires the 
use of density and internal energy floors. Because our calculation of 
the electron heating rate relies on quantifying the entropy changes 
in the underlying GRMHD solution, our predicted electron temper¬ 
atures also become unreliable when » pd. In the accretion disc 
simulations, this only affects the regions close to the pole where 
there is very little matter, not the evolution of the electrons in the 
bulk of the accretion disc or corona. We have specifically tested 


© 2015 RAS, MNRAS 000, 1-25 
















18 S. M. Ressler, A. Tchekhovskoy, E. Quataert, M. Chandra, C. F. Gammie 


several treatments of the internal energy and density floors which 
produce dramatically different results in the poles but are all con¬ 
sistent in the higher density regions for both the fluid variables and 
the electron temperature. 

Future applications of the methods developed in this paper 
will center on using our electron temperature calculations to pre¬ 
dict the emission from accreting black holes. In particular, we hope 
to produce more accurate images of the radio and IR emission of 
Sagittarius A* (and M87) that can be used to interpret the forthcom¬ 
ing spatially resolved observations by the Event Horizon Telescope 
(Doeleman et al. 2009) and Gravity (Gillessen et al. 2010). 
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APPENDIX A: TESTS OE ELECTRON CONDUCTION 

This Appendix outlines tests of our numerical implementation of 
electron conduction that demonstrate that our calculations are ro¬ 
bust and second-order accurate. The tests are taken directly from 
Chandra, Foucart & Gammie (2015), to which we refer the reader 
for more details. 


AI Conduction Along Eield Lines 


This test is simply to check whether the electrons conduct heat 
along field lines properly. The initial conditions are a 2D, periodic 
box of physical size 1x1 with uniform pressure and a small, density 
variation (and hence temperature variation) of the form; 

/ (v-0..'')^+(r-0.5)- \ 

P=Po(l-e )• (Al) 

The field lines are sinusoidal and given by 


S.t - Ba 

By = Bq sin(8;rx), 
derived from a scalar potential of 

1 


= Bo\y+ — cos(8;rx)|. 


(A2) 


(A3) 


For the conduction parameters, we choose^,, = 0.5/p and t = \ 
and run the simulation for IOt. 

Figure Al shows that the final state of the fluid is that of 
isothermal field lines, exactly as expected, with heat flux equili¬ 
brating the temperature along the magnetic field lines. This shows 
that our implementation of conduction properly limits the heat flux 
to be parallel to the magnetic field. 



0.2 0.4 0.6 0.8 


X 

Figure Al. Temperature profiles over-plotted with magnetic field lines 
in the 2D anisotropic conduction test from Chandra, Foucair & Gammie 
(2015), adapted for electron conduction (see §A1). The top panel is at the 
initial time while the bottom panel is at the end of the run (t = lOr). The 
field lines become isothermal, consistent with heat conduction only along 
the magnetic field. 


A2 Linear Modes Test 

This test checks whether our implementation of conduction gives 
the correct eigenmodes corresponding to Equation (19). Writing 
A = -a ± io), we initialise perturbations in an otherwise uniform 
box about the equilibrium solution with wave number k - 2 V2;r 
and run the simulation for one period: t = Inja). The analytic so¬ 
lution is that the perturbations, S, should obey S(t = 2;r/m) = S{t = 
0 )g- 2 ;ra/^ We choose S = 1/ V3x+ VI/ V3y and H = Inx + lny. We 
find that both (j) and converge at second order to the analytical 
solution as shown in Eigure A2. 

A3 ID Atmosphere in a Schwarzschild Metric 

This test checks whether our implementation of the electron con¬ 
duction gives the correct analytic result in a non-trivial space-time. 
In the Schwarzschild metric, the solution for a fluid in hydro-static 


equilibrium reduces to a system of two ordinary differential equa¬ 
tions, which can be solved for any given temperature profile (see 
Chandra, Eoucart & Gammie 2015 for details). Eor this test, we 
initialise the temperature and heat flux of the electrons to be this 
equilibrium solution for a purely radial field and see if the code 
can maintain it over a time of 100 r^/c in a computational domain 
of 1.4rg < r < 90rj. To compute the error, we again use the LI 
norm and find 2nd order convergence for both (p and as shown 
in Eigure A3. 

A4 Relativistic Bondi Accretion 

This test checks whether our implementation of the electron con¬ 
duction gives the correct analytic result in a fluid with u' 0, which 
activates terms that were not present in the ID atmosphere test. Eor 
the standard, spherically symmetric, steady-state Bondi solution for 
an accreting black hole (Hawley, Smarr & Wilson 1984), we can 
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Figure A2. LI Norm of errors in the 2D linear modes test after one period 
as computed from the eigenfrequencies given in equation (19). See §A2. 



Figure A4. LI norms of the error in the magnitude of the heat flux for the 
relativistic Bondi accretion test (§A4). 



Figure A3. LI norms of the error in both the heat flux and electron internal 
energy for the ID atmosphere test in the Schwarzschild metric (§A3). 


solve equation (16) by numerical integration if we assume that the 
heat flux does not back-react on the electron temperature. For this 
test, we set the initial condition of the fluid variables to be the Bondi 
solution and the initial conditions of (j) to he given by the solution 
to equation (16) with Dirichlet boundary conditions. We choose 
the sonic point to occur at = 20M and fix the outer boundary at a 
spherical radius of = 40M to have (j){r = 40M) = 0. The inner 
radius of the grid is inside the event horizon at r = 1.6M. The test 
is whether or not the code can maintain this state over a period of 
t = 200M. We find second order convergence of the heat flux to the 
analytical solution, as shown in Figure A4. 


APPENDIX B: DERIVATIONS 


BI Total Heating Rate 

This section derives the result quoted in equation (25). 

First, we introduce the variable kg, which is equivalent to Kg = 
PgP^^ at the beginning of the time step and at the n +1/2 “predictor” 
step, but which is evolved over a time step according to: 


dd = 0 . 


(Bl) 


We discretise equation (Bl) in a standard way (i.e. equation 21): 


At 

-,n+l/2 


rn«+l/2 


(B2) 


Ax 


where the square brackets indicate fluxes computed via the Rie- 
mann solver at cell interfaces and the generalisations to higher di¬ 
mensions is straightforward. Note that we have dropped the " in the 
n -I- 1 /2 and n terms because kg = Kg at the beginning of the time 
step and at the « + 1/2 step. We obtain the new value of entropy, 

at l„+i = t„+ At via solving equation (B2). We emphasise that 
kg is not the true entropy at but the entropy evolved according 
to equation (B1) [or its discretised equivalent equation B2] and thus 
does not include any heating. 

At the end of the time step (i.e. at f = f„+i), we compute the 
“true” value of the entropy due to the full GRMHD evolution, ac¬ 
cording to the definition of Kg'. 


K 


11+1 
g 



(B3) 


Unlike k"*^, which does not include any heating, Pg*' accounts for 
the heating as implied by the conservative evolution of the under¬ 
lying GRMFID scheme. The difference (Kg - kg)"*^ is related to the 
heating incurred during time step n, and we will use it below. 
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To compute the heating rate we evaluate the quantity: 

Q = pTifdfjSg - — — u'^dfjKg 

y -1 

= (B4) 

where the third equality holds because 


= d-gp^^)! ypg + pi^d^iKf 


and the first term vanishes due to conservation of mass. We evaluate 
eq. (B4) at the n + Ijl time step in a discretised form by centring 
the time derivatives at n + 1/2 but evaluating the prefactor at the 
n + 1 time step: 


-,,,+ 1/2 


1 pr -1 

7 - 1 yFs) 

(- ( yHSPKgU’f 


At 

,.xio+l /2 


(B5) 


lyF8PKgu%d -jyFs pKs'^'X 

Ax 


.vn„+l /2 


/ 1 \^+A /2 

Now, multiplying eq. (B2) by / yf-^j and adding the 

result to eq. (B5), we obtain equation (25) of the main text: 


-,n+l/2 


py ' {pu'(Kg-kg)]" 


At 


(B6) 


B2 Whistler Instability Limit on Conduction 


We assume that the electrons are relativistic with 0,, = kTcIm^c- > 
1. If the electrons relax to thermal equilibrium with a scattering rate 
Vj, relativistic kinetic theory implies that the electron viscosity rp, 
and thermal diffusivity satisfy (Anderson & Witting 1974) 

A A 

rje-Qe— Xe- 1-6— (B7) 

Velocity space instabilities set an upper limit on the electron ther¬ 
mal conductivity in a turbulent plasma. Physically, as the magnetic 
field in the accretion disc fluctuates in time, this generates pres¬ 
sure anisotropy, which is resisted by velocity space instabilities 
that isotropise the distribution function and thus limit the magni¬ 
tude of the thermal diffusivity. Chandra et al. (2015) show that the 
theory of relativistic anisotropic viscosity implies that the pressure 
anisotropy and scattering rate are related by 


AP, 

V, —— =: In 

p 




(B8) 


where AP,, = P± - Py and we have neglected some general rela¬ 
tivistic terms for simplicity. 

Electrons satisfy limits on pressure anisotropy of 


AP, 1.3 AP, 0.25 


(B9) 


The second term on the right hand side of equation (B9) is a fit 
to the whistler instability threshold for relativistically hot electrons 
(based on numerical solutions of the dispersion relation derived in 
Gladd 1983). The coefficient in the numerator technically depends 


weakly on 0,, varying from 0.125 for non-relativistic electrons to 
0.25 for 0, =: 10 (Lynn 2014). Note that the slope of the /3, term 
for the whistler instability in equation (B9) is a fit for,8, - 0.1 - 30. 
Gary & Wang (1996) and Sharma et al. (2007) found a somewhat 
shallower slope oc in non-relativistic calculations but this is 
not a good fit over the large dynamic range of ,0, considered here. 
The first limit in equation (B9) is the electron firehose instability 
which is an electron-scale resonant analogue of the fluid firehose 
instability (Gary & Nishimura 2003). This limit is based on non- 
relativistic calculations and should to be extended to the relativistic 
limit in future work. However, based on our whistler calculations 
this is unlikely to be a significant effect. 

Sharma et al. (2007) found that the typical pressure anisotropy 
satisfied AP/P > 0 in simulations that explicitly evolved a pressure 
tensor. Physically, this sign of the pressure anisotropy corresponds 
to outward angular momentum transport. Assuming that the RHS 
of equation (B8) is ~ G, the whistler instability limit in equation 
(B9) thus implies,^, ~ crj(r/r^)^^^(4j0,)“° ®. This is not a significant 
constraint on the conductivity relative to the saturated value ixe ~ 
cvg), for ,0, < 1, which can occur either in the corona/outflow or 
because P, <K Tp. However, this estimate does suggest that the 
electron conductivity may be modest in the bulk of the disc at 
lOtg ify0, » 1. 

Equation (B9) can be implemented by calculating AP,/P, us¬ 
ing equation (B8) given an assumed(and using equation (B7) to 
relate v, and;^^,). If equation (B9) is violated, v, should be increased 
and;^^, decreased such that equation (B9) is satisfied. Alternatively, 
an even simpler first approximation would be to simply limit 

A. :S crg{rlrgf\4lt,r°-^ = (BIO) 

motivated by the estimate in the preceding paragraph for the 
whistler instability. This is the limit we have used in the accretion 
disc simulations in §5 of the main text. 


B3 Electron Conduction Numerical Stability 

Non-relativistically, an explicit implementation of thermal con¬ 
duction is stable only if the time step. At, satisfies the condition 
At < Ax^ lx, where Ax is the grid spacing in 1-dimension and,y is 
the thermal diffusivity. The relativistic theory outlined in §2.3, how¬ 
ever, where the heat flux tj) is evolved according to equation (16), 
differs from the non-relativistic case in that it is no longer diffu¬ 
sive. This alters the criterion for stability to be a condition on the 
relaxation time, t, given by equation (36), which we derive here. 

To check the numerical stability of our conduction theory we 
assume that we are in Minkowski space in the rest frame of the 
fluid, and further simplify our analysis to one dimension in which 
b = i. 

Under these assumptions, a Von Neumann stability analysis 
on equations (16) and (5) leads to a quadratic equation for the am¬ 
plification factor, G, with the following solutions: 


G = 1 - C[l -cos(fe)] 


1 At 


2 T 


Xe-t 


--(l±,/l-4(y,-l)|^sin(/:) 


(BID 


with the condition for stability being that |G| < 1. Here, as before, C 
denotes the Courant factor. To analyse equation (Bll), we consider 
two cases: 1) when the square root term is real, and 2) when the 
square root term is imaginary. 

When the square root term is real, the condition for stability 


© 2015 RAS, MNRAS 000, 1-25 













22 S. M. Ressler, A. Tchekhovskoy, E. Quataert, M. Chandra, C. F. Gammie 


becomes: 


T > l^t 


2 - C [1 - cos(k)] - (j, - 1)^ sin (kf 


(B12) 


(2-C[l -cos(fe)])2 

The right hand side is a maximum for k = n modes, which gives, 
simply: 

At 


T > 


2 ( 1-0 


(B13) 


The more interesting case is when the square root term in 
Equation (Bll) is imaginary, where the criterion for stability be¬ 
comes: 


T ^ At 


(r. - 1)^ sin {kf -H C [1 - cos(A:)] - 1 


C(2 - C[1 - cos(yt)]) [1 - cos(fc)] 
which, defining K = {je — \)AtXelAA , has a maximum at 


cos(k) 


c- V 4 ff(i - 0-0 

a^-2K{l-C) 


if 


At > 


Ax^ 1 - 4C(1 - O + VI - 4C(2C - 1) 


Ge - lk<. 

= At 


8 ( 1-0 


(B14) 


(B15) 


(B16) 


and a maximum at fc = ;r otherwise. So if At < our criterion 
becomes: 


T > At 


2C- 1 


4C(1 -O 

Finally, if At > At„i,, then we have 

T > 

2A:(0-4C(1-0) 




(B17) 


Af X 


4KC(l -O(V41^(l - 0-0 - 2C) -H 2C^ 
(40(1 -O + 4KC{\ -O -O) V401 - 0-0 


4KC(\ -O(V401 - 0-0 - 2C) + 2C^ 


(B18) 


— 'fmax,2‘ 

The general behaviour of Equation (B18) is complicated, but the 
result is roughly consistent with 

I At 

’U 

for most reasonable choices of the Courant factor. This is the result 
quoted in equation (36) of the main text. 

To summarise, our scheme is stable when: 


T > 


At 


2 ( 1-0 

At 


I Tjiflx, 1 


, 'tmax,2 


: At < Atr, 


: At > Atr, 


L 2 (l -0 

for At„i,, Tmax.i, and as defined in equations (B16), (B17), and 
(B18), respectively. 


B4 Electron Heating in a ID Shock 

Formally, for an ideal shock in a zero-viscosity fluid there is no 
unique path in (P,p) space that connects the pre and post-shock 
values given by the Rankine-Hugoniot conditions, meaning that the 
dissipation per unit volume, JpTds, is not a well-defined quantity. 


However, by introducing any non-zero viscosity, the degeneracy 
is broken and there exists a unique path in (P,p) space and hence 
a well-defined dissipation. To see this, we take the ID Rankine- 
Hugoniot relations for a static shock, given some prescription for 
the viscous stress, t = 4l3p^ -v (pis the dynamic viscosity coeffi¬ 
cient, and can be an arbitrary function of plasma parameters). 


m - pv 


P 

E 


■ pv^ + P + T 
1 


(B20) 


, 7 

—pv H- -Pv + TV, 

2 7 — 1 


where m, p, and E are constants representing the mass, momentum, 
and energy flux across the shock. Absent t, we could combine these 
three equations in several different ways to get a relationship of the 
form P = P(p). With non-zero viscosity, however, there is only one 
unique way to do this, namely, by taking pv- E and solving for m, 
which gives: 


llirP E 

P(P) = (7- 1) r- P+-P 

\l p in 


or, in terms of k = Pp 


I \ / 1 ^ P ^ 


(B21) 


(B22) 


We assume that the electrons receive a constant fraction of the total 
heat: 


pTelPdf,Se - fepTgU'^d^Sg 


o' o' 

-ifdgiKe = fe - rl/duKg 


(B23) 


1 


“7-1 


or in quasi-conservative form (using the mass continuity equation 
and assuming a flat space metric): 


- {pu'^Kr) = fg— -- (pU^Kg) . (B24) 

ax" ^ ■' r - r ax" 


The final electron entropy is given by integrating this equation from 
the initial to the final density, which, for a ID shock reduces to 


/ 


(niKe) - fe 


Pf 

7.-1 f 

— J' 


p^ ^'m-—dp. 
dp 


(B25) 


giving: 

= n' 


K = u‘A — 
Pi 


fe I in j + 1 _ ^ 7_ + T ~ 1 
7 — 1 \ 2p/ "Xe + 1 7e lit 7e — 1 


(B26) 


/. Ipf 
7- 1 \Pi 


m J + 1 j Epi 7 — 1 

2 p, ye + 1 ^ 7 . m je-I 


where pf is determined from the Rankine-Hugoniot conditions. For 
a strong shock with Mach number » 1, this simplifies to 


: mVi 


fe 


Ji-l 


7 -t- 1 

7- 1 


i- 21 ) + i + Z 


(B27) 


where u, is the pre-shock fluid velocity in the shock’s rest frame. 
Dividing by Ug = 2mVj ( 7 ^ - 1^ yields equation 41. 
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APPENDIX C: ELECTRON HEATING IN A VISCOUS 
SHOCK 


In this appendix we show that by introducing an explicit bulk vis¬ 
cosity to the non-relativistic hydrodynamic equations, our electron 
heating calculation outlined in §3.3 give an electron internal en¬ 
ergy that converges to the analytic result derived in Appendix B4 
for electron heating at a shock. 

We treat viscosity by explicitly adding the ID viscous energy 
and momentum fluxes to the ideal MHD fluxes for a constant kine¬ 
matic viscosity, v: 


4v du 

E,vise — 

3 ax 
4v du 

Fp.visc = - 


(Cl) 


Note that these are non-relativistic fluxes which are formally incon¬ 
sistent with the relativistic code in which they are used. However, 
our goal here is simply to show that with a resolved shock struc¬ 
ture the electron heating calculation converges to the correct an¬ 
swer. The non-relativistic limit is fine for this purpose. The fluxes 
in equation (Cl) smooth out discontinuities to a continuous profile 
of finite width, determined by v and the velocity scale. The solution 
for the profile of a viscous shock, now defined as a smooth transi¬ 
tion from an initial to final state as opposed to a discontinuity, can 
be computed analytically for a constant kinematic viscosity, v. In 
the shock frame, taking jc ^ -oo as the initial state, this solution 
takes the form: 


v(x) 


(^+^l)q.exp 

3(x-xo)Pi/ 1\ 

4v \ M) 

exp 

3(x-xo)v,/j 1\ 

4v \ Mj 

+ (7 + 1) 


(C2) 


where M is the pre-shock Mach number, w, is the pre-shock speed 
al X -oo, and xq is a constant determining the location of the 
shock. For a pre-shock density p,, the density profile is obtained 
from the mass conservation equation: pjvjv(x), which determines 
the pressure profile from equation (B21). Similarly, the profile for 
the internal energy of the electrons in terms of p(x} is given by 
equation (B26) with the substitution p/ ^ pix). 

For our numerical test, we do not use the standard Noh test 
as outlined in §4.1.2 due to the problems noted by the original pa¬ 
per (Noh 1987). For any numerical scheme that gives the shock a 
finite width, the formation of the shock from the converging flow 
undershoots the density at the center of the grid by a finite amount 
that does not disappear at higher resolution. Given this difficulty, 
our numerical test is instead to set the initial and boundary con¬ 
ditions of both the fluid and electron variables equal to the ana¬ 
lytic solution for a stationary shock (e.g., equation C2) and evolve 
for a dynamical time of L/u,, where L is the grid size. We choose 
7 = 5/3, Vi = lO^^c, M ~ 49, and v = O.Olw/L. Figure Cl shows 
both the density profile and the ratio of the electron internal en¬ 
ergy to the total internal energy for both jg = 4/3 and = 5/3 
electrons at the end of the run as compared to the analytic solution 
(equation B26). We find good agreement with the analytic solution 
and second order convergence (Figure C2) up to the resolution at 
which relativistic errors in the analytic solution become important 
{6uglug ~ (u/c)^ ~ lO^"*). 



X 



X 


Figure Cl. High Mach number (~ 49), stationary, viscous shock results 
for an electron heating fraction fg = 0.5 at a resolution of 2000 cells. Top: 
fluid density. Bottom: electron internal energy relative to total fluid internal 
energy. Both the jg = 4/3 and yg = 5/3 electrons display good agreement 
with the analytic solution, converging at 2nd order (see Figure C2). This is 
in contrast to the formulation without explicit viscosity used in §4.1.2, in 
which the shock structure is always just a few grid points. An accurate cal¬ 
culation of the shock heating requires a well-resolved shock structure (i.e., 
a shock with a finite width), which is provided by adding explicit bulk vis¬ 
cosity to the fluid equations. Given that the error incurred by our numerical 
scheme without explicit viscosity (~ 3%) is acceptable for our purposes, we 
do not use explicit viscosity in our calculations. 
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Figure C2. Convergence results for the electron internal energy in a steady- 
state, ID, high Mach number, viscous shock as compai'ed to the analytic so¬ 
lution (see Appendix C). Both the je = 5/3 and je = 4/3 electrons converge 
at 2nd order, as opposed to the non-viscous shock of §4. 1.2 where only the 
je = 5/3 electrons converged to the analytic solution. Second order con¬ 
vergence is achieved in this test problem because the shock profile is well- 
resolved and continuous. This shows that our method correctly captures the 
dissipation in strong shocks when the shock profile can be resolved. At the 
highest resolution, relativistic corrections to the (non-relativistic) analytic 
solution become important so the error no longer converges at second or¬ 
der. 
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APPENDIX D: TORUS INITIAL CONDITIONS 

In this appendix we describe in more detail the initial configuration of the torus in our simulations of an accreting black hole. In all expressions 
that follow we measure radii in units of the gravitational radius Vg = GMjc?- (or equivalently set G = M = c = 1). 

Fishbone & Moncrief (1976) derived an equilibrium solution (their equation 3.6) of the general relativistic hydrodynamic equations in 
the Kerr metric in terms of the relativistic enthalpy, h = (p + Pg + Ug)lp, and the constant angular momentum per unit mass, I = Ug,u' . We use 
their equation 3.6 exactly as presented when r > rj„ and when the right-hand side is positive, otherwise we set p = P = 0. Additionally, we 
assume an adiabatic equation of state, P - kqp^, for some choice of kq, and fix / such that the density maximum occurs at 


l=< 


2,Cl '^rjjiax "1" 

d 1 

' maxj 

1 1 

(«^-2aVw + r^ax)] 

V 

2a V'^max + '"max 

— s y* 

-''max 


+ '■max - 2rmax) yjI +2arJJ^ -3Ir 


(Dl) 


yt-max + '■^ax “ 3''max 


where a is the dimensionless spin parameter of the black hole. This expression for / is equivalent to the Keplerian value at r = rn, 
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